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ABSTRACT 

We compare the results of various cosmological gas-dynamical codes used to simulate 
the formation of a galaxy in the ACDM structure formation paradigm. The various 
runs (thirteen in total) differ in their numerical hydrodynamical treatment (SPH, 
moving-mesh and AMR) but share the same initial conditions and adopt in each case 
their latest published model of gas cooling, star formation and feedback. Despite the 
common halo assembly history, we find large code-to-code variations in the stellar 
mass, size, morphology and gas content of the galaxy at z = 0, due mainly to the 
different implementations of star formation and feedback. Compared with observation, 
most codes tend to produce an overly massive galaxy, smaller and less gas-rich than 
typical spirals, with a massive bulge and a declining rotation curve. A stellar disk is 
discernible in most simulations, although its prominence varies widely from code to 
code. There is a well-defined trend between the effects of feedback and the severity 
of the disagreement with observed spirals. In general, models that are more effective 
at limiting the baryonic mass of the galaxy come closer to matching observed galaxy 
scaling laws, but often to the detriment of the disk component. Although numerical 
convergence is not particularly good for any of the codes, our conclusions hold at two 
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different numerical resolutions. Some differences can also be traced to the different 
numerical techniques; for example, more gas seems able to cool and become available 
for star formation in grid-based codes than in SPH. However, this effect is small 
compared to the variations induced by different feedback prescriptions. We conclude 
that state-of-the-art simulations cannot yet uniquely predict the properties of the 
baryonic component of a galaxy, even when the assembly history of its host halo is 
fully specified. Developing feedback algorithms that can effectively regulate the mass 
of a galaxy without hindering the formation of high-angular momentum stellar disks 
remains a challenge. 

Key words: cosmology: theory - methods: numerical - galaxies: formation - galaxies: 
evolution 



1 INTRODUCTION 

Numerical simulations play a central role in studies of cos- 
mic structure formation. CoUisionless N-body simulations 
have now become the main theoretical tool to predict the 
non-linear evolution of dark-matter dominated structures 
once initial conditions are specified. Their high accuracy and 
huge dynamic range have allowed a detailed comparison of 
their outcome with observations of the large-scale structure 
of the universe. The impressive agreement between these ob- 
servations and the predictions of the A Cold Dark Matter 
(ACDM) model has helped to e stablish it as the cu rrent 
paradigm of structure formation l|Springel et al.|[2006l ). 

Simulating the evolution of the visible Universe is 
much more complex, because it requires understanding 
the many astrophysical processes which drive the evolu- 
tion of the baryonic component under the gravitational 
influence of the dark matter. Numerical hydrodynam- 
ics in cosmological simulations has traditionally used ei- 
ther the Lagran gian Smoot h ed Particle Hydrodynamics 



technique (SPH ; (LucvI 1 19771: Icingold fc MonaghanI 1 19771 : 



iMonaghanl 119921 : iKatz et al T ll996l : ISpringell l2010bl ) or Eu- 
lerian grid-based solvers sometimes aided by Adaptive 



1992 
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Mesh Refinement (AMR) techniques (|Cen fc Qstriked 
Brva n fc Norman" 199 51: lKravtsovlll999l : iFrvxell et ali 
Tevssicr 2002; Quilis 20041'). 

Both approaches have advantages and disadvantages. It 
is widely appreciated that SPH is not able to capture shocks 
with high accuracy and that in certain situations fluid in- 
stabilities can be supp ressed , at least for standar d imj)le- 
mentations of SPH (iAgertz et al.ii2007i : iCreasev et al.ii201H ) . 
On the other hand, mesh-based codes are not Galilean in- 
variant and may in some cases generate entro py spuriously 
through artificial mixing (IWadslev et al.ll2008h . As a result, 
even for some simple non-radiative problems, Lagrangian 
and Eulerian codes do not converge to the same solution 
(e.g., lOkamoto eralll2003l : IAgertz et allbOOTl : iTasker et ahl 



12008 



Mitchell et al.ll2009l l. Novel techniques, such a s the La- 



yangia n, moving- mesh AREPO code introduced bv lSpringell 
2010ah . hold the promise of improving this state of affairs, 
but their application is still in its infancy. 

An even more uncertain ingredient of galaxy formation 
simulations are the descriptions of poorly understood physi- 
cal processes such as star formation and feedback. The huge 
dynamic range between the super-Megaparsec scales needed 
to follow the hierarchical growth of a galaxy and the sub- 
parsec scales that govern the transformation of its gas into 



stars implies that direct simulation of all relevant physical 
processes is still out of reach of even the most powerful com- 
puters and best available algorithms. Star formation and 
feedback are therefore introduced in cosmological simula- 
tions as "sub-grid" parameterized prescriptions of limited 
physical content and lacking numerical rigor. 

These difficulties have hampered the progress of sim- 
ulations of galaxy formation within the ACDM paradigm, 
but a few results of general applicability have nevertheless 
emerged. For example, the formation of realistic disk galax- 
ies in dark matter halos formed hierarchically, as expected 
in ACDM, was re cognized as a challenge ev en in early simu- 
lation s (see, e.g.. iNavarro et al. 1995: Na varro fc Steinmet j 
I1997I I. Simulated galaxies suffered from "overcooling" and 
from a dearth of angular momentum due to the trans- 
fer of angular momentum from the baryonic component 
to the dark matter during the ma ny merger events that 



characterize hierarchi c al ass e mblv jNavarro fc Benz 



ISommer-Larsen et al.1 fl999l : 1 avarro fc Steinmet j I2OOGI ) 



I991I: 



Feedback, as a general heating mechanism that can pre- 
vent overcooling and regulate the assembly of a galaxy 
whilst avoiding catastrophic angular momentum losses, 
emerged as a crucial ingredient of any successful galaxy for- 
matio n simulation (IWeil et al.l 19981: Navarro fc Steinmetj 



1997"; 'Sommcr-Lar sen et al.lll999l : IScannapieco et al.ll2008l : 
Zavala ot al. 2008). 

Considerable progress has been made since this time: 
recent simulations have shown that the angular mo- 
mentum problem can indeed be alleviated when feed- 
back from evolvin g stars, in parti c ular supernovae (SNe), 
is inc luded (e.g., Okamoto et al.l 20051: Governato et al.l 



20071: IScannapieco et all I2OO8I: ICeyerino fc KlvpinI I 2OOG 



Scannapi eco et al.l2009l:|joung et al.|2009l : IColm et al.i20ld : 
Sales et al.l l2010l : [stinson et al.ll201Gl ). Some authors have 



also investigated alternative feedback mechanisms, such 
as energy liberated during the formation of supermas- 
black holes as well as th a t carried by cosmic 



Jubolgas ot al. 


2008; Booth & Schayc 


2009: Fabian et all 


2OIOI: IWadepuhl fc Springellboill). 



Despite progress, difficulties remain. For example, sim- 
ulations tend to allow far too many baryons to accrete into 
galaxies to be consistent w ith the observed stellar mass func- 
tion of galaxies (see, e.g.. lGuo et al.ir201Gl ). In addition, al- 
though stellar disks do form, they are often too concen- 
trated, with stee ply-declining rotatio n curves at odds wit h 
observation (e.g.. lA.badi et aLll2003al lbl: IStinson et al.ll2010l ). 
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The difficulties in simulating disk galaxies highlighted 
above are compounded by the limited guidance afforded 
by analytic and semi-analytic models of galaxy forma- 
tion. In such modeling, the properties of disk galaxies 
are typical ly envis ioned to r eflect t hose of their surround- 



ing halos (iFall fc Efstathio ul Il980l: lKa,uffmann et aT 



Dalcanton et al.l Il997i: Mo et all Il998l : iBower et al 



199a 



200e 



De Lucia fc Blaizoj 2007 ): for example, halos with high net 
spin and quiet recent merger histories are commonly as- 
sumed to be likely sites of disk galaxy formation. How- 
ever, there are growing in dications that these ass ump- 
tions might b e too s implistic. IScannapieco et ahl l|2009l ) and 
IStinson et al.l (|2010l ). for example, find no clear relation be- 
tween the presence of disks and the spin parameter of the 
halo: disk s form in halos with low an d high spin parameters. 
Moreover, IScannapieco et al] (|2009l ) report that disk forma- 
tion is not assured by a quiescent assembly history, since 
they can be destroyed not only by major or intermediate- 
mass mergers, but also by secular processes such as a mis- 
alignment between the gaseous and stellar disks. 

These difficulties have led to little consensus on what 
determines the morphology of a galaxy; what the main feed- 
back mechanisms are; and what role they play on different 
mass scales and at different times. Indeed, there is even de- 
bate about whether the difficulties in reproducing realistic 



merical resolution fe.E.. iGovernato et al.l 


20041), inappropri- 


ate modeling of the relevant physics (e.g.. 


Maver et al.ll2008l: 



model (e.g., Sommer-Larsen fc Dolgovll200ll ). 



Progress in this unsettled field requires at least a careful 
evaluation of the different numerical techniques, resolution, 
and choice of sub-grid physics adopted by various authors. 
Most groups do carry out and publish resolution tests and 
convergence studies of their own numerical setup. However, 
because of the complexity of the problem, the lack of telling 
test Ccises with known solution, and the absence of clear 
theoretical predictions for individual systems, each group 
chooses to tune the relevant numerical parameters accord- 
ing to different priorities and/or prejudice, and there has so 
far been little effort invested in comparing the results of dif- 
ferent techniques and codes. Would they give similar results 
if they followed the formation of a galaxy in the same dark 
matter halo? 

The main goal of the Aquila Project is to address this 
question by comparing the predictions of different codes us- 
ing common initial conditions and a homogeneous set of 
analysis tools. Rather than focusing on whether individual 
codes perform better or worse than others, we contrast their 
predictions for the stellar mass, angular momentum content, 
star formation rates and galaxy size, with observation. 

This paper is organized as follows. § [2] describes the 
initial conditions and the simulation setup. § [3] compares 
the galaxy morphology, star formation history, size, angular 
momentum, and gas content of the 13 simulated galaxies, 
together with a brief discussion of the effects of numerical 
resolution on the results. Finally, § |4] summarizes our main 
findings and lists our main conclusions. 



2 THE SIMULATIONS 
2.1 The Codes 

The Aquila Project consists of 13 different gas-dynamical 
simulations of the formation of a galaxy in a ACDM halo of 
similar mass to that of the Milky Way. Nine different codes 
were used for the project (two codes were run three times 
each with different sub-grid physics modules). The various 
codes differ in their hydrodynamical technique (SPH, AMR, 
moving-mesh): seven codes use the SPH technique, six of 
which are based on GADGET (hereafter G3 for short, see 

I ! ' ' ' 

[Springel 2005) and one on GASOLINE (hereafter referred to 
as GAS: Wadslev et al.ll2o"oi ). 

The G3-based codes share the same numerical grav- 
ity /hydrodynamical treatment but differ in their cool- 
ing/s tar formation/feedback modules. G3 refers to the stan- 
dard ISpringel fc Hernguistl l\200^ implem entation; G3-CS 
refers to the code presented in Scannapieco et al.l 
200(|; G3-T O to that developed by lOkamoto et al.l 



I200a bold ): G3-GIMIC i s described in Grain et al . ., 
G3 -MM is introduced bylMurante et al.l (|2010l ): and G3-CK 




bv lKobavashi et al. I(|2007l). Of the two codes that do not use 
SP H one is the r amses AMR code (hereafter R, for short, 
see lTevssieijl200^). and the other is the moving-mesh code 
AREPO (|Springelll20"l0al ). 

Each code was run by the group responsible for its de- 
velopment, adopting (independently from the choices made 
by other Aquila participants) their latest published model 
of cooling, star formation, and feedback. These differ from 
code to code. Regarding radiative cooling, for example, some 
codes assume primordial abundances to compute cooling 
rates; others use metal-dependent cooling rate tables; and in 
some cases cooling is implemented on an element-by-element 
basis. Star formation also varies from code to code, although 
in nearly all cases the efficiency of transformation of gas 
into stars is set by attempti ng to reproduce the Kennicutt- 
Schmidt empirical relation (lKennicuttlll998l ) in simulations 
of isolated disks (see also Fig. lAip . 

The numerical treatment of feedback also varies from 
code to code. In most cases, thermal feedback is used, where 
supernova energy is injected into the interstellar medium 
(ISM) as thermal energy. The dispersal of the input energy 
is usually delayed artificially, in order to promote the pres- 
surization of the ISM, the onset of winds, and the effective 
regulation of subsequent star formation. A few codes adopt 
kinetic feedback, where kinetic energy is dumped directly 
into the gas. In the G3-TO code, the outflowing gas is tem- 
porarily decoupled hydrodynamically from the rest of the 
ISM to ensure that the specified mass loading (wind mass 
per unit mass of stars formed) and velocity are not modified 
by viscous drag until gas escapes the star forming region. 

Two further simulations were run with the G3 code: 
one (G3-BH) that included, in addition to SN, the feedback 
energy associated with the assembly of supermassive black 
holes and a third (G3-CR) where another form of feedback, 
that associated with energy deposition by cosmic rays, was 
also included. 

Two further RAMSES runs are also part of our series, 
one (R-LSFE) where the star formation timescale is much 
longer than the fiducial choice, delaying substantially the 
transformation of gas into stars, and another (R-AGN) where 
the feedback energy was augmented by the contribution of a 
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Figure 1. Maps of the dark matter distribution in the region surrounding the Aquila halo at 2 = 0. The dark matter is projected in 
cubic volumes of side length as indicated in each panel. Pixel brightness corresponds to the dark matter density using a logarithmic scale. 



putative AGN associated with a central supermassive black 
hole. 

The main characteristics of each code are summarized 
in Table[T] Appendix|X]presents a more detailed description 
of each of the codes and their numerical choices. 

2.2 Initial Conditions 

All simulations share the same initial conditions (ICs), a 
zoomed-in resimulation of one of the halos of the Aquar- 
ius P roject (halo "Aq-C", in the notation of ISpringel et al.l 
l2008f ). Th e ICs were Rene r ated u sing Fourier methods as de- 
scribed in ISpringel et al.l (|2008l ) , modified to include a gas 
component. The displacement field is first calculated on a 
set of grids and then interpolated onto the nodes of the un- 
perturbed particle positions, chosen from a glass-like config- 
uration. For SPH codes these particles represent the matter 



distribution, while for the AMR initial conditions they rep- 
resent the dark matter only. 

The displacement field is used to perturb the particle 
positions and to assign them velocities consistent with the 
growing mode of the density fluctuations. For AMR runs the 
density and velocity flelds of the gas are needed on a set of 
meshes. These quantities were calculated in the same way 
as described above except that the displacement and density 
fields were interpolated onto the vertices of a pre-defined set 
of nested hierarchical grids tailored to requirements of the 
AMR code. 

For SPH runs, each high-resolution dark matter parti- 
cle is split into two to create one dark matter and one gas 
particle, with relative masses given by the assumed value 
of the universal baryon abundance parameter fib (Table [2]). 
The positions of the new particles are such that their center 
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Table 1. Summary of code characteristics and sub-grid physics 



Code 


T3 icifViT-iCiri f 


Tvno 


UV background 
(^Uv) (spectrum) 






G3 (gadgets) 


[1] 


SPH 


6 


[10] 


primordial [13] 


SN (thermal) 


G3-BH 


[1] 


SPH 


6 


[10] 


primordial [13] 


SN (thermal), BH 


G3-CR 


[1] 


SPH 


6 


[10] 


primordial [13] 


SN (thermal), BH, CR 


G3-CS 


[2] 


SPH 


6 


[10] 


metal-dependent [14] 


SN (thermal) 


G3-TO 


[3] 


SPH 


9 


[11] 


element-by-element [15] 


SN (thermal-f kinetic) 


G3-GIMIC 


[4] 


SPH 


9 


[11] 


element-by-element [15] 


SN (kinetic) 


G3-MM 


[5] 


SPH 


6 


fiol 


primordial [13] 


SN (thermal) 


G3-CK 


[6] 


SPH 


6 


[10] 


metal-dependent [14] 


SN (thermal) 


GAS (gasoline) 


[7] 


SPH 


10 


[12] 


metal-dependent [16] 


SN (thermal) 


R (ramses) 


[8] 


AMR 


12 


[10] 


metal-dependent [14] 


SN (thermal) 


R-LSFE 


[8] 


AMR 


12 


[10] 


metal-dependent [14] 


SN (thermal) 


R-AGN 


[8] 


AMR 


12 


[10] 


metal-dependent [14] 


SN (thermal), BH 


AREPO 


[9] 


Moving Mesh 


6 


[10] 


primordial [13] 


SN (thermal) 



NOTE : [1] Springel et"all ll2008l') : [2] IScannapieco et al ll20q5h : IScannapieco et alj ll2006l 'l: f3]IOkamoto et alj ll201(t: f4llCrain et alJ 
1I2OO9I'): f5][Mu rante et al.' |'2010'); [6] 'Kobav ashi et al.1 h oO?)-. f7]'Stinson et al.' 1*2006'): [81 'Tevssiei^ l'2002l'): iRasera fc Tevssied (120061 ): 
iDubois fc Tws sier (2008); [9] Springel (2010£^'l: flOl lHaar dt & Madau (1996); [11] Haardt & Madau (2001);J12] Haa rdt fc Madau (200 5, 
private communication); [131 IKatz et alj (1 199^ ; [141 [Sutherland fc Dopitai (Il993l 'l: [151 IWiersma et al.l ((2009^); [16l lShen et al.l 1 I2OIOI ). 



Table 2. Code parameters for simulations at level-5 resolution (level-6 
parameters are given between parentheses) 



Code 


/b 


muM 


^gas 


Softening 






(Qb/an) 


[lO^Mo] 


[10«Mo] 


£|=° [kpc] 


Zfix 


G3 












G3-BH 












G3-CR 


0.16 


2.2 


0.4 


0.7 





G3-CS 




(17) 


(3.3) 


(1.4) 


(0) 


G3-CK 












Arepo 












G3-TO 


0.18 


2.1 


0.5 


0.5 


3 


G3-GIMIC 




(17) 


(3.7) 


(1) 


(3) 


G3-MM 


0.16 


2.2 


0.4 


0.7 


2 






(17) 


(3.3) 


(1.4) 


(2) 


GAS 


0.18 


2.1 


0.5 


0.46 


8 






(17) 


(3.7) 


(0.9) 


(8) 


R 


0.16 


1.4 


0.2 


0.26 


9 


R-LSFE 




(11) 


(1.8) 


(0.5) 


(9) 



R-AGN 



NOTE: /b: baryon fraction; rriDM- mass of dark matter particles in the 
high resolution region; m,gas: initial mass of gas particles; ^g~^- gravi- 
tational softening at 2 = 0; za-^: redshift after which the gravitational 
softening is kept fixed in physical coordinates. The softening is fixed 
in comoving coordinates at 2 > zgx (see Appendix [Bt . 
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• #' 




10 kpo 



A * 




Hk.„„ = 1,40110" Ms M.,,u.r = a.!)!)-*!;!" 




E-LSPE 






M^u.r = S,l9>!10" Ma 




H^u.r = 3.S3110" Ms 



projected mass densiLy [log(MQ / kpc^)] 



6.50 7.00 7.50 B.OO 8.50 9.00 9.50 10.00 10.50 



Figure 2. Face-on and edge-on maps of projected stellar mass density. The face-on projection is along the direction of the angular 
momentum vector of galaxy stars. The face-on and edge-on maps are 30 X 30 kpc^ and 30 X 12 kpc^ , respectively. The size of each pixel 
is 58.6 pc on a side and its color is drawn from a logarithmic color map of the surface stellar mass density. The total stellar mass within 
the galaxy radius (rg^i =0.1 r2oo ~ 25 kpc) is listed in the legend of each panel. 



of mass position and velocity are identical to those of the 
original particle. 

The selected halo, Aq-C, has a present- day mass simi- 



lar to the Milky Way 1.6 x 10^^ Mr) (e.g. iDehnen et all 



20061: iLi fc White! l2008l : ISmith et al.ll2007l : IXue et al.ll200a 



Watkins et all |2010 '~ and has a relatively quiet formation 
history. It is also mildly isolated at 2: = 0, with no neighbor- 
ing halo more massive than half its mass within a radius of 
1 Mpc. Maps of the dark matter distribution in boxes 
of various sizes are shown in Fig. [T] 



2.3 Cosmology 

We assume a ACDM cosmology with the following parame- 
ters: ilm = 0.25, f^A = 0.75, (Tg = 0.9, Us = 1, and a Hubble 
constant oi Ho = 100 /i kms"^ Mpc"^ = 73kms"^Mpc"\ 
These parameters are consistent with the WMAP 1- and 
5-year results at the 3a level and are identical to the pa- 
rameter s used for the Millenn i um and Millennium-II simu- 
lations l|Springel et al.ll2005bl : iBovlan-Kolchin et al.ll2009l '). 
The value of f2b used in each simulation is given in Tabled 

The Millennium-II is, in fact, a resimulation of the cos- 
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A 




-1.6 -1.0 -0.5 0.0 0.5 1.0 1.5 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 -1.5 -1.0-0.5 0.0 0.5 1.0 1.5 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 




-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 



R-LSFE 





-1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 -1.5-1.0 -0.5 0.0 0.5 1.0 1.5 -1.5-1.0-0.5 0.0 0.5 1.0 1.5 -1.5-1.0-0.5 0.0 0.5 1.0 1.5 



Figure 3. Distribution of stellar circularities, e = jz/jc, for the different models. The circularity parameter is the z-component of the 
specific angular momentum of a star particle, j^, expressed in units of the circular orbit value, jc, at that radius. Stars with e 1 typically 
belong to a rotationally-supported disk component. Thick and thin lines correspond to level-5 and level-6 resolution runs, respectively. 



mological volume from which the Aquarius hales were origi- 
nally selectecJB- Aq-C is thus prese nt both in t h is sim ulation 
and in the semi-analytic model of IGuo et al.l |20l3) which 
was tuned to fit the luminosity, stellar mass, size and gas 
content functions measured fo r galaxies in the Sloa n Digital 
Sky Survey (SDSS). Similarly, [Cooper et al.l (|2010l ) used the 
GALFORM code to model star formation in all six Aquarius 
halofl using parameters very similar to those o nBower et alj 
l|2006i) . which reproduce local galaxy luminosity functions. 
The properties found for the central galaxy of Aq-C in these 
two models thus give an indication of what direct simula- 
tions should produce if implementation on a cosmological 
volume is to reproduce observed galaxy abundances. 



2.4 The Runs 

All 13 simulations were run at two different numerical reso- 
lutions: level 5 and level 6, respectively, following the naming 
convention of the Aquarius project (lower numbers indicate 



^ In the convention of the Aquarius Project (lower numbers in- 
dicate higher resolution), the Millennium-II simulation has a res- 
olution intermediate between levels 5 and 6 (the two resolutions 
used in our set of simulations). 

^ In this case, the level-4 resolution simulations were used. 



higher resolution). Table [5] gives the dark matter and (ini- 
tial) gas particle masses for the two resolutions. The gravi- 
tational softening is kept fixed in comoving coordinates until 
redshift zax, after which its value is fixed in physical coordi- 
nates. The simulations vary in their choice of zax and there- 
fore the gravitational softening has slightly different values 
at 2 = 0, listed as e|=° in Table [5] (see also Fig. IBT]|. 



2.5 Analysis 

We describe here some of the conventions and definitions 
used in the analysis of the simulated galaxies. The center 
of the galaxy is defined to coincide with the position of the 
baryonic particle with minimum gravitational potential. The 
virial radius, r2oo, is the radius of a sphere, centered on the 
galaxy, with mean density equal to 200 pcrit , where pcrit = 
3H^/8tvG is the critical density for closure. We use the term 
halo to refer to all the mass within r2oo and galaxy to the 
baryonic component within a radius rgai = 0.1 x r2oo from 
the center. 

Where a distinction is drawn between "hot" and "cold" 
gas, we adopt a temperature threshold of lO^K to sepa- 
rate the two phases. Some codes (e.g. G3, G3-BH, G3-CR, 
G3-GIMIC, G3-TO, AREPO) adopt an "effective" equation of 
state to describe the ISM and to circumvent numerical in- 
stabilities in poorly-resolved regions. This may cause some 
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fluid elements to have nominal temperatures in excess of 
10^ K, but still be star- forming. In order to prevent assign- 
ing this gas to the hot phase, we automatically assign all 
star-forming gas particles to the "cold" phase. 

As we shall see below, different runs yield simulated 
galaxies of widely-varying baryonic mass and angular mo- 
mentum. In particular, not only the specific angular momen- 
tum changes between simulations (as expected, given the 
wide range in galaxy mass spanned by the various runs), but 
also its orientation, as we discuss in Appendix [C] Because 
of this, for orientation-dependent diagnostics, we rotate each 
simulated galaxy to a new coordinate system where the an- 
gular momentum vector of its stellar component coincides 
with the 2 direction. 



3 RESULTS 

We present here results concerning the stellar mass, mor- 
phology, size, star formation history, and angular momen- 
tum content of the simulated galaxy. Unless otherwise spec- 
ified, all results correspond to z — and to the level-5 resolu- 
tion runs. Numerical convergence between level-5 and level-6 
runs is discussed in 



3.1 Galaxy morphology 

Fig.[2]shows face-on and edge-on maps of the projected stel- 
lar mass density for the 13 runs. Labels in each panel list 
the simulation name, as given in Table[T] as well as the total 
stellar mass of the galaxy (i.e., within rgai). 

These figures illustrate the complex morphology of the 
simulated galaxies; bars, bulges, and extended disks are 
present, but their relative prominence varies widely from 
run to run. The galaxy stellar mass also shows large scatter, 
spanning about a decade from the least (G3-TO) to the most 
massive (R), respectively. 

A quantitative measure of the importance of a 
rotationally-supported component is provided by the distri- 
bution of stellar circularities, e, defined as the ratio between 
the 2-component of the specific angular momentum of a star 
and that of a circular orbit at the same radius r: 



Jz ^ Jz 

ic(r) rK(r)' 



(1) 



where Vc{r) = GM{< r)/r is the circular velocity. Stars 
belonging to a disk are expected to have e ~ 1, whereas stars 
belonging to a non-rotating spheroidal component should 
hav e an e-distribution roughly symmetric around zero (see, 
e.g.. lAbadi et al.ll2003bl : IScannapieco et al.ll2009l '). 

We show the circularity distribution of all 13 runs in 
Fig. [3l Thick and thin lines correspond to the level-5 and 
level-6 resolution simulations, respectively. The diversity in 
morphology seen in Fig. [2] is clearly refiected in the distri- 
bution of circularities. Thin disks that appear prominently 
in the images show up as well-defined peaks in the circu- 
larity distribution at e ~ 1, a distinction that sharpens at 
higher numerical resolution. In some cases, notably G3, G3- 
MM, G3-CK, and AREPO, the galaxy is noticeably fiattened 
and clearly rotating, but lacks a prominent thin disk. 

The importance of a thin disk may be crudely estimated 
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Figure 4. Median formation time of stars in the galaxy at z = 0, 
expressed in terms of the expansion factor, 050% = 1/(1 + -250%): 
as a function of the fraction of stars with circularity exceeding 
0.8. 



by the fraction of stars with e > 0.8, /(e > O.sfl Only in 
four simulated galaxies do more than ~ 40% of stars sat- 
isfy this condition, two SPH-based and two AMR-based: R, 
R-LSFE, G3-GIMIC, and GAS. The most extreme case, R- 
LSFE, provides a clue to this behaviour. In this simulation 
feedback is inefficient and star formation is deliberately de- 
layed, allowing gas to accrete into the galaxy and settle into 
a centrifugally-supported structure before turning into stars. 

Indeed, any mechanism that hinders the early transfor- 
mation of gas into stars without curtailing gas accretion later 
on is expected to promote t he formation of a disk (see, e.g., 
iNavarro fc Steinmet j [l997l ) . As a result, the galaxies with 
most prominent disks a re also the ones with the youngest 
stars (|Agertz et al.|[201ll ). This is shown in Fig.[3], where we 
plot /(e > 0.8) versus the median formation time of all stars 
in the galaxy (expressed in terms of the expansion factor, 
0,50%)- A. clear correlation emerges, with disks increasing in 
prevalence in galaxies that make their stars later. On the 
other hand, galaxies that make their stars early tend to be 
spheroid-dominated. 

An interesting outlier to this trend is G3-MM, which 
forms stars as late as R but has a small fraction of stars in 
a disk. Further investigation shows that the G3-MM galaxy 
did harbour a disk, but it was severely impacted by a col- 
lision with a massive satellite in recent times. The satellite 
is present in other runs, but it has not yet collided with the 
main galaxy in the majority of cases. This is due to the fact 
that even small differences in the early evolution get am- 
plified with time and can lead to large discrepancies in the 
orbital phase of satellites later on. To the extent that this 
can infiuence the morphology of the central galaxy, a certain 
degree of stochasticity in the morphological evolution of a 
simulated galaxy seems unavoidable. 



Note, however, that these fractions often compare poorly with 
photometric estimates of the disk-to-total ratios dAbadi et alj 



.2003a : .Scannapieco et alllioiol) . 
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Another interesting result to note is that neither G3 
nor AREPO form prominent thin disks. These two runs share 
the same sub-grid physics, but use very difTerent numerical 
hydrodynamical techniques, which suggests that the mor- 
phology of the simulated galaxies is indeed rooted mainly 
in how gas gets accreted and transformed into stars and in 
the merger hist ory of the pa r ticula r halo. As discussed re- 
cently by, e.g., iTorrev et al.l l|201ll 'l. the numerical scheme 
does make a difference when considering the detailed prop- 
erties of simulated disks, but it does not seem to be the 
main reason why the G3 and AREPO runs lack disks in this 
halo. Rather, the failure of feedback to regulate effectively 
the onset of early star formation and to allow for late gas 
accretion seems the most likely culprit. Support for this in- 
terpretation is provided by G3-BH and G3-CR which, despite 
the increased feedback, fail to prevent most stars from form- 
ing quite early. As a result, none of these models allows a 
sizeable thin disk to develop. 

Models with more efficient feedback schemes, such as 
those where feedback regulates more effectively early star 
formation through galactic winds (e.g., G3-CS, G3-TO, G3- 
GIMIG, gas), yield galaxies with two well-defined compo- 
nents: an old, non-rotating spheroid surrounded by a young 
rotationally-supported disk. Still, even in this case the disk 
component is subdominant in terms of total stellar mass, 
with /(e > 0.8) around ~ 30-40%. 

Finally, it is worth noting that the morphology of a 
galaxy is often dissimilar at the two resolutions attempted 
here. In general, more prominent disks form at higher reso- 
lution but in some cases this trend is reversed. We further 
discuss resolution effects in Section [TSl 



It seems one could argue that a successful feedback mech- 
anism must selectively remove l ow-angular momentum ma- 
terial from the galaxy (see , e.g.. iNavarro fc Steinmet j|l997l : 
Ivan den Bosch et al.ll20o3 : iBrook et al.ll201ll )! 



3.2 Rotation curves 

As discussed above, regulating star formation without hin- 
dering the formation of a stellar disk is a challenging task 
for any feedback implementation. One might argue that a 
solution might simply be to delay star formation, such as 
in R-LSFE, but this comes at the expense of unrealistic disk 
properties. A simple and convincing diagnostic is the "rota- 
tion curve" of the disk which, for simplicity, we represent by 
the circular velocity profile of the galaxy, Vc{r). 

This is shown in Fig. [S] where we group in four panels 
the results of the 13 level-5 Aquila runs, and compare them 
with the circular velocity curve of the dark matter-only Aq- 
C run (dark solid line) and, for referen ce, with the rotati on 
curve of the Milky Way as compiled bv lSofue etakl (|2009l fl. 

This figure makes clear that the "best disks" in terms of 
morphology (i.e., R-LSFE, R, GAS, and G3-GIMIC) aU have 
steeply declining rotation curves, at odds with the flat ro- 
tation curves characteristic of normal spirals. The extreme 
R-LSFE model again illustrates the problem: here feedback 
is inefficient at removing baryons, allowing large amounts of 
gas to collect in a central disk before being turned into stars. 
A large fraction of these baryons have relatively low angular 
momentum, however, leading to the formation of a disk that 
is unrealistically concentrated and with a declining rotation 
curve. (A similar consideration applies to G3 and AREPO.) 



" Note tliat lSofue et al.l | |2009| '| assume a Galacto-centric position 
and velocity of tlie sun of 8 kpc and 200 km s~^, respectively. 



Support for this view comes from inspection of the rota- 
tion curves of galaxies where galactic winds play a more sub- 
stantive role, especially at high-redshift, when low-angular 
momentum baryons are preferentially accreted: G3-CS and 
G3-TO show nearly flat rotation curves. A similar result is 
found for G3-BH, G3-CR and R-AGN, but in these cases the 
"success" must be qualified by noting that none of these 
galaxies have a clearly-defined stellar disk: the fiat Vc-curves 
are just a refiection of the low baryonic mass of the galaxy 
that results from adopting these extremely effective feedback 
models. 



3.3 Stellar Mass 

The stellar mass of a galaxy is determined by the combined 
effects of radiative cooling, the rate at which cold gas is 
transformed into stars, and the ability of feedback to reg- 
ulate the supply of star-forming gas. Fig. [B] shows how the 
various implementations affect the stellar mass of the central 
galaxy, Mstdiar. This figure shows Mstdiar as a function of 
^^200 (2), the virial mass of the main progenitor from z = 2 
to jz = 0. (The symbols correspond to values at z = 0.) 

To guide the interpretation, we show with a dashed 
curve the total baryonic mass within the virial radius cor- 
responding to the universal baryon fraction, (f2b/J^m).A/2oo, 
which sets a hard upper limit to the stellar mass of the 
central galaxy. The shaded region surrounding the dotted 
curve corresponds to the stellar masses predicted, at 2: = 0, 
by requiring agreement between the halo mass and galaxy 
stellar mass func tions through simple abundance matching 
(|Guo et al.ll20ld ). We also show the prediction for Aq-C of 
the semi-ana lytic models GAL FORM (Cooper et al.ji201Q ) and 
L-GALAXIES (|Guo et aLlbOlll l. 

The most striking feature of Fig. [6] is the large code-to- 
code scatter in the stellar mass of the galaxy, which varies 
between ~ 4 x 10^° Mq and ~ 3 x lO" Mq at z = (Ta- 
ble[3|. The three largest stellar masses are obtained with the 
mesh-based codes, R, R-LSFE and AREPO, and correspond to 
assembling nearly all available baryons in the central galaxy. 
This illustrates the weak efficiency of the feedback imple- 
mentations chosen for these codes, aided by the fact that, 
at comparable resolution, cooling efficiency is enhanced in 
mesh-based codes relative to SPH jVogelsberger et alll201ll : 
ISiiacki et af]|201ll : lKeres et al.ll201lf r 

Indeed, G3 forms only about half as many stars as 
AREPO, despite sharing exactly the same sub-grid physics. 
The galaxy formed by GAS also has a large stellar mass, but 
this may be due to the fact that this code has a more ef- 
ficient cooling function through the addition of metal-line 
cooling. Although the numerical technique may effect some 
changes in the stellar mass, these are small compared with 
the variations introduced by the feedback implementation. 
This may be seen by noting that, when including AGN feed- 
back in RAMSES, the stellar mass decreases by a factor of ~ 5 
and the disk component is largely erased. 

Note that the large variations in the stellar mass of the 
galaxy formed in different runs imply that the dark halo will 
respond by contracting differently in each case, as we discuss 
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Figure 5. Circular velocity curves of all galaxies in the level-5 runs, for the resolution level-5 runs. The four panels group the results 
according to numerical technique. The top-left panel corresponds to various feedback choices of the standard GADGET code; the top-right 
and bottom-left correspond to other, independent star formation/feedback modules developed for GADGET, as well as the SPH-based 
GASOLINE code. The bottom-right panel groups the results of the AMR code RAMSES and the moving-mesh code AREPO. Thick and thin 
lines correspond to level-5 and level-6 resolution runs, respectively. The solid circles indicate, for the level-5 simulations, the position of 
the stellar half-mass radius of each simulated galaxy. The thick black line shows the circular velocity of the dark-matter-only simulation 
of the same halo (Aq-C). For reference, the region shaded in light grey is bounded by the peak and viri al velocities of the A quarius halo. 
Dark grey points with error bars are observed data for the Milky Way's rotation curve, as compiled bv lSofue et al. I ||2009D . 



in Section [3.41 Note also that the differences in stellar mass 
are dominated by differences prior to 2; = 2; in fact, in some 
simulations the stellar mass at 2; = 2 is already above the 
z = stellar mass-halo mass relation. 

It is also important to note that feedback must be 
roughly as effective as that of R-AGN in order to obtain stel- 
lar masses consistent (within the error) with the abundance- 
matching predictions. Indeed, the only other codes to match 
this constraint, and thus fall within the shaded area of Fig. [6] 
are G3-BH and G3-TO; of these only the latter forms a galaxy 
with a discernible disk (see Fig. [2}. All other models give 
stellar masses well in excess of the abundance-matching con- 
straint, a shortcomi ng of most publ i shed galaxy format ion 
simulations to date (|Guo et al.ll2010l : ISawala et al.ll201ll ). 

It is also worth noting that the abundance-matching 
models allow for substantial scatter in the M2oo-Afsteiiar re- 
lation. Indeed, the more sophisticated treatments of the L- 
GALAXIES and GALFORM semianalytic codes indicate that 
Aq-C might form a galaxy more massive than expected on 
average for a halo of that mass (see open and filled starred 
symbols in Fig. [6|. l-GALAXIES, in particular, suggests that 



Aq-C might be a 2a outlier from the relation, which would 
alleviate, but not resolve, the disagreement between the re- 
sults of R, R-LSFE, AREPO, and GAS and the model pre- 
dictions. GALFORM, on the other hand, predicts that Aq-C 
should be about la above the mean abundance-matching 
relation. 

Taken altogether, these results illustrate the basic chal- 
lenge faced by disk galaxy formation models: feedback must 
be efficient enough either to prevent the accretion, or to fa- 
cilitate the removal, of most baryons, whilst at the same 
time allowing enough high-angular material to accrete and 
form an extended stellar disk. 



3.4 Tully-Fisher relation 

The stellar mass and circular velocity of disk galaxies 
are strongly linked by the Tully-Fisher relation, and it is 
therefore instructive to compare the properties of simu- 
lated galaxies with those of observed disks. T his is done in 
Fig. [71 where we compare data compiled by iDutton et al.l 
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Figure 6. The stellar mass of the central galaxy as a function of the virial mass of the surrounding halo. Curves of different color track 
the evolution of the galaxy in each simul ation betwe e n z = 2 and ^ = 0. The dotted line indicates the stellar mass expected at 2 = 
from the abundance-matching analysis of I Guo et alj | |201C| ): the shaded region corresponds to a 0.2 dex uncertainty. The dashed line 
indicates the mass of all baryons wi thin the virial radiu s, (Hb/^m) M200 . The filled and open star symbols indicate the predictions of 
the semi-analytic models GALFORM l|Cooper et al.ll20ld 'l and L-GALAXIES ijGuo et al.ll201lf l for halo Aq-C, respectively. The dot-dashed 
curves show the evolution since z = 2 according to these two models. 



(l201lf) from IPizagno et all |2003), IVerheiienI l|200ll l. and 
ICourteau et 3.1.1 ^ 2007 ) with the 13 simulated galaxies. 

Because the rotation curves of simulated galaxies are 
not flat (see Fig. [SJ, we use velocities estimated at the stel- 
lar half-mass radius in order to be as consistent as possi- 
ble with the rotation speeds estimated obse rvationally from 
spatially-resolved rotation curves (see, e.g.. ICourteau et al.l 
[2007). The symbols connected by a solid line show the con- 
tribution of the dark matter to the circular velocity at the 
same radius. A dotted line shows the same, but for the dark 
matter-only Aq-C halo; the difference between solid and dot- 
ted curves indicates the degree of "contraction" of the dark 
halo. 

There is a clear discrepancy between the observed TuUy- 
Fisher relation and simulated galaxies, which tend to have 
substantially larger velocities at given Afstdiar- The disagree- 
ment worsens for large stellar masses, emphasizing again 
the fact that too many baryons are able to cool and form 
stars in these systems. Interestingly, at low stellar mass 
simulated galaxies approach the observed relation but still 
have, on average, higher rotation speeds than typical disks. 
This suggests that, although these galaxies may have stellar 
masses consistent with abundance-matching considerations 
(see 13. 3|) . they must differ from typical spirals in other re- 
spects, such as an excessive concentration of the dark matter 
or luminous component. 

The dark matter contribution to the circular velocity 
(connected symbols in Fig. [7]) lies well below the average 
rotation speed expected from the TuUy-Fisher relation. This 



suggests that the concentration of dark matter is not the 
origin of the disagreement; there should in principle be no 
problem matching the observed relation provided that the 
luminous component of the galaxy is extended enough. The 
offset from the observed TuUy-Fisher relation thus suggests 
that simulated galaxies are more concentrated than normal 
spirals, resulting in disks that rotate too fast for their stellar 
mass. We analyze the size of simulated galaxies in §[3]6l after 
examining the importance of the gaseous component of the 
galaxy next. 



3.5 Gaseous component 

Fig. [8] shows /gas, the fraction of the baryonic mass of sim- 
ulated galaxies in form of gas at 2 = 0, as a function of 
the R-band absolute magnitude, and compares the m with 



data for star forming galaxie s from Schombert et al.l ((200 ll ) , 
iBell fc de Jon3 (|2000|) and iHavnes et all (|l999l ). Magni- 



tudes for the simulations w ere calculated using the dust-free 
iBruzual fc Chariot] (|2003 l) population synthesis models, for 
a Chabrier Initial Mass Function (IMF) and solar metallic- 
ity. 

Most simulated galaxies have gas fractions below 10%, 
which puts them at odds with observations of nearby spirals. 
As for the stellar mass, note the large code-to-code scatter 
in /gas, which varies from about 1% for G3 to nearly 30% 
for R-LSFE. As expected, galaxies with larger gas fractions 
are predominantly those with morphologies that include a 
well-defined stellar disk, presumably because of ongoing star 
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Figure 7. The TuUy-Fisher relation. Tiie circular velocity at the stellar half-mass radius of each simulat ed galaxy is plotted a s a function 
of stel lar mass for all 13 level-5 r uns. Small black dots correspond to data for nearby spirals taken from lPizaeno et alj l(20o3) , l^rheiie^ 
1I2OOII ) and lCourteau et alj ||2007| '|. The symbols connected by a solid line show the contribution of the dark matter to the circular velocity 
at _Rii stars- Those connected by the dotted line show the circular velocity of the dark matter-only halo (Aq-C) at the same radii. 



formation. The converse, however, is not always true: G3-CS 
and R have low gas fractions a± z = but prominent disks. 

More surprisingly perhaps, the gas fraction seems to 
correlate only weakly with the present-day star formation 
rate (which we discuss more thoroughly in § 13. 7[) . For ex- 
ample, R-AGN has the third largest gas fraction and by far 
the lowest star formation rate at 2; = 0. The same applies 
to G3-TO, which, despite its large /gas, forms stars at rates 
well below what would be expected for an average spiral (see 
Fig.IH]). 

Overall we see no obvious dependence of the gas fraction 
on the numerical method: of the four galaxies with highest 
/gas, two are SPH-based (G3-TO and G3-MM) and two are 
AMR-based (R-AGN and R-LSFE). However, we note that 
AREPO has a much higher gas fraction (and stellar mass) 
than G3, despite sharing the same sub-grid physics. This 
supports the conclusion that standard SPH-based methods 
may underestimate the total amount of gas that cools and 
becomes available for star formation, especially when feed- 
back is as weak as implemented in the G3 and AREP O runs 
(see also lAgertz et aLll2007l : IVogelsberger et al.llioTll ') . 

The interpretation of these results is not straightfor- 
ward. The gas content of a galaxy is constantly evolving; 
supplied by accretion, depleted by star formation, and re- 
moved by feedback-driven winds. This leads to large fluctu- 
ations in the instantaneous gas fraction and star formation 
rate of a galaxy, which may be exacerbated by chance events 
such as satellite accretion. 



3.6 Galaxy size 

The left-hand panel of Fig. [5] compares the stellar half-mass 
radii of simulated galaxies with the size of observed galax- 
ies of similar stellar masfl, as well as with the predictions 
of the GALFORM and l-GALAXIES semi-analytic models and 
the Milky Way, for reference. The observed sizes correspond 
to Petrosian half-light radii in the r-band rather than stellar 
half-mass radii so the comparison should only be taken as in- 
dicative. Red/blue dots correspond to galaxies redder/bluer 
than (g-r) ^ 0.59-1-0.052 logio(A/steiiar/M0) - 10.0 and are 
meant to outline roughly the location of early- and late-type 
galaxies in this plot. 

As anticipated in the previous subsection, most galaxies 
are too compact to be consistent with typical spiral galax- 
ies of comparable stellar mass. In general, the more massive 
the simulated galaxy, the smaller its size, a trend that runs 
contrary to observation. In particular, the most massive sim- 
ulated galaxies (R, R-LSFE, AREPO) are even smaller than 
most early- type galaxies in the nearby Universe, highlight- 
ing again the shortcomings of simulations where cooling and 
star formation proceed largely unimpeded by feedback. 

Simulations where feeedback is more effective at curtail- 
ing the formation of stars give rise to galaxies with sizes in 
better agreement with observation. For example, the half- 
mass radii of R-AGN and G3-BH reach _Rh ~ 5 kpc for 



5 Data taken from the SDSS MPA-JHU DR7 release for nearby 
(z < 0.1) galaxies; http://www.mpa-garching.mpg.de/SDSS/DR771 
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Figure 8. Gas mass fraction, /gas = Mgas/(Mgas + Mstollar), of the galaxy versus R-band absolute magnitude. Magnitudes have been 
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Figure 9. Left: Stellar half-mass radius as a function of the stellar mass of the galaxy. Blue and red dots show the Petrosian r-band half 
light radius for a sample of nearby (z < 0.1) SDSS galaxies taken from the MPA-JHU DR7 release. The sample is split into "blue cloud" 
and "red sequence" galaxies depending on their colors, according to the condition (g — r) = 0.59-1- 0.052 login(^Btellar/-^w) ~ 10-0 (only 5 
per c ent of randomly se lected data point s are shown) . We also show the predictions of the semi-analytic models GALFORM l lCooper et alj 
l201Ch and l-GALAXIES iIGuo et al.ll20lil ) for Aq-C, and the approximate location of the Milky Way in this plot. Right: The projected 
half-mass ra dius of cold (T < 10 ^ X) g as in the galaxy as a funct ion of stel l ar ma ss. Grey circles indicate the half-mass radii of HI disks 
compiled bv lDutton et aL I 1I2OIII') from lSwaters et al.l l l 19991) and I Verheiieiil 1I2OOII '). and the open star symbol indicates the prediction of 
L-GALAXIES. The prediction of GALFORM is not shown here, since it predicts a present-day gas mass close to zero. 
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Figure 10. Median formation time of stars in the galaxy (ex- 
pressed in terms of tlie expansion factor (050%) as a function 
of total stellar mass at 2 = 0. We also show the pred iction of 
the semi-a nalytic models GALFORM ijCooper et aljBoiOh and L- 
GALAXIES jGuo et al.ll201ll 'l for Aq-C. 



~ 7 X 10^" Mq, at the lower end of the distribution of spiral 
sizes. However, as discussed in i|3.1l neither of these galaxies 
has a disk-like morphology. 

The simulated galaxy with lowest Mstoiiar and a well- 
defined disk is G3-TO, but, as seen from Fig.O it has a half- 
mass radius of less than 2 kpc, well below what would be 
expected for a spiral of that mass. This is because most of the 
stellar mass in G3-TO is in the form of a highly-concentrated 
spheroid rather than in an extended disk. Therefore, even in 
this case, feedback has apparently allowed too many low 
angular momentum baryons into the galaxy. 

These results support our earlier conclusion: feedback 
must not only limit how many baryons settle into the galaxy, 
but must also selectively allow high angular momentum ma- 
terial to be accreted and retained in order to form a realistic 
disk. 

The right-hand panel of Fig. [5] shows the projected 
half-mass radius of cold gas in the simulated galaxies (at 
2 = 0) as a function of stellar mass and compares them 
with HI observation s co mpiled by Dutton et al.l (|201ll ') from 
ISwaters et"all (|l999l ) and lVerheiienI (|200ll '). Despite the large 
code-to-code variation, the simulated gaseous disks are sys- 
tematically more extended than the stellar component, in 
agreement with observation. They are also in better agree- 
ment overall with the typical size of HI disks, a result that 
suggests that material accreted relatively recently (and thus 
still in gaseous form) has, on average, enough angular mo- 
mentum to form disks of realistic size. If feedback were able 
to favor the accretion and retention of this late-accreting, 
high-angular momentum gas, then simulated galaxies would 
have a much better chance of forming realistic disks. 



3.7 Star formation history 

A recurring theme of our discussion so far has been the need 
to prevent the assembly of an overly massive galaxy with- 



out, at the same time, preventing high-angular momentum, 
late-accreting material from reaching the galaxy and form- 
ing a disk. This requires feedback to act especially strongly 
at high redshift, when a large fraction of baryons first be- 
come cold and dense enough to start forming stars. None 
of the Aquila runs seems able to meet these requirements 
satisfactorily, as shown by the star formation history of the 
simulated galaxies. 

Fig. 1101 shows the stellar mass of the galaxy versus the 
median formation time of the stars (expressed in terms of 
expansion factor, a — 1/(1 + z)). Note the strong correla- 
tion between the two, which indicates that the codes best 
able to limit the growth of the mass of the galaxy do so at 
the expense of curtailing the incorporation of late-accreting 
material. Indeed, the three galaxies with the lowest Mstciiar 
(G3-TO, R-AGN, smaU G3-BH) form half of their stars in the 
first Gyr or so of evolution, i.e., by z ~ 4. It is not surprising 
then that two of them lack a discernible disk, and that the 
disk in G3-TO is overwhelmed by a massive, dense spheroid 
composed mainly of old stars. 

Further details on the star formation history of indi- 
vidual galaxies are presented in Fig. 1111 where we show, in 
cumulative and differential form, the distribution of stellar 
ages of the stars in the galaxy at 2 = 0. Note that these are 
not star formation rates for the main progenitor, since stars 
may (and some of them, indeed, do) form in different pro- 
genitors before being accreted into the galaxy. Nevertheless, 
the data in Fig. [TT] show clearly that few codes are able to 
prevent the early burst of star formation activity that ac- 
companies the collapse of the first massive progenitors of the 
galaxy. Only G3-TO, G3-MM, R-AGN and R-LSFE are suc- 
cessful at keeping this peak "rate" at less than ~ 100 Mq 
yr~^ at z ~ 4-5, but even they see a precipitous decline in 
star formation afterwards. (The exception is R-LSFE, but 
this is achieved by artificially delaying star formation, see 

§[n]) 

Fig. [TT] also shows that the morphological appearance 
of simulated galaxies is very weakly correlated with star for- 
mation history. There are examples of galaxies with lots of 
recent star formation that have well-defined disks (e.g., R) 
and examples which do not (e.g., AREPO), as well as cases 
of galaxies with very little recent star formation that have 
well-defined disks (e.g., G3-CS, G3-T0) and cases which do 
not (R, G3-BH). 

Aside from these general considerations, the details of 
the star formation history of each galaxy reflects the partic- 
ular suh;;gnd_2hysicsimplem^ chosen for each code 
(see ISchave et all I2OI0I ) . For example, the feedback scheme 
of G3-GIMIC (where SN-driven winds are assumed to be 
launched with fixed speed) imply that it is effective at early 
times, when the mass of the progenitors is small, but be- 
comes ineffective when the main halo reaches ~ 10^^ Mq . 
This curtails early star formation but allows the stellar 
mass to increase substantially at low redshift. On the other 
hand, feedback in G3-TO is effective at all redshifts, since its 
strength scales with the potential well of the galaxy. In G3- 
CS star formation is not effectively regulated at very early 
times because SN energy feedback is not injected into the 
ISM instantaneously but rather after a time delay which 
depends on the local properties of the cold and hot neigh- 
bouring particles. 

These choices imply that at any given time there is sub- 
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Time [Gyr] Time [Gyr] Time [Gyr] Time [Gyr] 

1 3 5 7 10 13 1 3 5 7 10 13 1 3 5 7 10 13 1 3 5 7 10 13 




Figure 11. Distribution of stellar formation times (expressed in terms of the expansion factor, a). Left and right panels show the same 
data, in differential and cumulative form, respectively. The simulations are grouped according to numerical technique, as in Fig. \E\ The 
squares and circles indicate the time of formation of the first 10% and 50% of the stars. 

stantial scatter in the star forming properties of simulated 
galaxies, compounded by the fact that there is a certain de- 
gree of stochasticity in the rate at which a galaxy accretes 
mass 13. ip . For example, in the case of AREPO, the infall 
of a satellite at 2; ~ 0.7 disrupts the gas disk and leads to 
a significant increase in the star forming activity at 2 = 0. 
In the GAS run, a large burst of star formation also occurs 
near z = 0, but in this case it seems associated with en- 
hanced gas accretion facilitated by the infall of a satellite. 
Such effects are at least partially responsible for the large 
code-to-code scatter in the star formation rate of the galaxy 
at the present time. This is shown in Fig. 1121 where we com- 
pare the present-day star formation rate (averaged over the 
past 0.5 Gyr to smooth out short-term fluctuations) with the 
stellar mass of simulated galaxies and constrast them with 
observations of local SDSS galaxies. The observational sam- 
ple corresponds to nearby {z < 0.1) SDSS galaxies selected 
from the MPA-JHU DR7 catalofl The SFR of simulated 
galax;ies varies from a low of ~ 0.02 M© yr~^ (R-AGN) to 
nearly 20 yr^^ (gas), spanning the whole range covered 
by observed galaxies, from "red and dead" spheroids to ac- 
tively star-forming gas-rich disks. The large scatter leads us 
to conclude that caution must be exercised when analyzing 
the instantaneous star formation rates of simulated galaxies, 
since these depend sensitively on the details of accretion and 
of the implemented sub-grid physis. 



3.8 Numerical Convergence 

Convergence is a necessary (but not sufficient) criterion to 
assess the robustness of numerical simulation results. We dis- 
cuss here the effects of resolution by comparing the results 

^ |http : //www . mpa-garchlng ■ mpg . de/S DSS/DR7/ 1 
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Figure 12. Present-day star formation rate (averaged over the 
last 0.5 Gyr to smooth out short-term fluctuations) as a func- 
tion of stellar mass. Blue and red dots correspond to a sample 
of nearby {z < 0.1) SDSS galaxies selected from the MPA-JHU 
DR7 catalog (only 5 percent of randomly selected data points are 
shown) . The sample is split into "red sequence" and "blue cloud" 
galaxies as described in Fig. [9] We al so show th e pred iction of 
the semi-analytic model L-GALAXIES of lCuo et al.l 1I2OIII) for Aq- 
C and the a pproximate location of the Milky Wav llOliver et all 
l20ld iLeitnc r & Kravtsov 2011). Note that the GALFORM semi- 
analytic model is not shown here since it predicts a present-day 
star formation rate close to zero. 



of level-5 and level-6 simulations for each code (see also Ta- 
ble [3]). As a quick reminder, at level 5 the parent halo of the 
Aquila Project, Aq-C, has ~ 700, 000 dark matter particles 
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Figure 13. Numerical convergence. We show tiie fractional variation in various galaxy properties between the level-5 and level-6 
resolution runs of each code: AQ/Q = {Qg — Q^j/Q^. Results are shown for the stellar mass (A/gtciiar) a-nd half-mass radius (iJh,stars); 
for the peak circular velocity (Vmax) and corresponding radius (rmax); for the gas mass (Afgas) and half-mass radius (i?h,gas); for the gas 
fraction (/gas) and present-day star formation rate (SFR); and for the median star formation time ((150%) and the fraction of stars with 
circularities larger than 0.8, /(e > 0.8)). Arrows indicate that results lie outside the plotted range. The actual values of the quantities 
used for the plot are listed in Table [S] 
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within the virial radius; this number is reduced to ~ 90, 000 
at level 6. 

Also, when considering convergence one must note that 
resolution-dependent behaviour is inevitable to some degree, 
since the Jeans length and Jeans mass of star forming gas are 
not well resolved either at level 6 or at level 5. For instance, 
in several codes stars are allowed to form only above a den- 
sity threshold of nth ~ 0.1 cm~^ and at temperatures of 
order ~ 10^ K, which correspond to a Jeans length of ~ 1.5 
kpc, only slightly larger than the gravitational softening at 
level 6. 

Although all codes start from initial conditions of the 
same mass and spatial resolutions, each code then adopts 
its own refinement strategy, which may result in significant 
effective runtime resolutions. For example, some codes of 
similar type choose different gravitational softening lengths, 
and differ as well in their choice of when and whether to 
keep it fixed in comoving or physical coordinates (Fig. IBl| . 
This makes it even harder to disentangle numerical resolu- 
tion effects from code type effects. The enhanced cooling in 
grid-based codes discussed in Sec. 13.31 for example, seems 
driven by differences inherent t o the hydrodynamical tr eat- 
ment, and not by resolution (see lVogelsberger et al .I2OIII . for 
details), but the distinction is less clear for other quantities. 

We quantify the effects of resolution on a given quantity, 
Q, by the fractional variation, measured at z = 0, between 
the level-5 and level-6 runs, 

AQ/Q = (Qs - Q5)/Q5. (2) 

Note that with this definition a quantity that increases with 
increasing resolution will have a negative vale of AQ/Q. 

We summarize the results in Fig. 1131 where we show 
the variations: (i) in global galaxy properties such as peak 
circular velocity, Vmax and the radius at which it is achieved, 
fmax; (ii) in properties of the stellar component, such as the 
total mass, Mstciiar and the half-mass radius, -Rh, stars; (iii) 
in properties of the gas component, such as the mass, A/gas, 
and half- mass radius, -Rh.gas; (iv) in properties related to 
star formation, such as the gas fraction and the present- 
day star formation rate; and (v) in the details of the galaxy 
assembly /morphology, such as the median formation time of 
stars at z — (expressed in terms of the expansion factor, 
150%) and the fraction of stars kinematically associated with 
a rotationally-supported disk (/(e > 0.8), see 13. ip . When 
differences exceed 100% (i.e., \AQ/Q\ > 1), we use arrows 
to indicate if such deviations occur along the x and/or y 
coordinates of the plot. The actual values of all quantities 
plotted in Fig. [T3] are listed in Table O 

Aside from the fact that numerical convergence is not 
particularly good for any of the codes. Fig. [13] illustrates a 
few interesting points. The first concerns the properties of 
the stellar component (left two upper panels). In particu- 
lar, the total mass in stars and their median age seem to 
be the most reliable results, suggesting that the star for- 
mation/feedback scheme chosen for each code is reasonably 
independent of resolution. Most codes reproduce the total 
stellar mass and a^o% to better than 20-30%. Interestingly, 
the peak circular velocity is one of the properties least af- 
fected by resolution (see also Fig. [5]). This is encouraging, 
since it implies that diagnostics such as the observed TuUy- 
Fisher relation may be used to assess the success of a par- 
ticular model. 



Convergence deteriorates when considering the detailed 
properties of the stellar component, such as the fraction of 
stars in high-circularity orbits: /(e > 0.8). Although half 
the codes give results that converge to better than ~ 30% 
the variations can be much larger for some codes. Inter- 
estingly, increasing the resolution does not always leads to 
better-defined disks: for example, the fraction of "disk" stars 
increases by a large fraction in the case of GAS but actually 
decreases for AREPO. Finally, there is some indication that 
the most extreme feedback models (i.e., those that result in 
the lowest Mstoiiar; e.g., R-AGN and G3-BH) are the most 
vulnerable to resolution-induced changes. 

Even larger variations are seen for the properties of the 
gas component: only five of the thirteen simulations show 
variations in Mgas and -Rh.gas smaller than 50%, and three 
simulations have differences in A/gas larger than 100 percent. 
Similar results are found for the gas fractions, and, conse- 
quently, for the present-day average star formation rates. 
In general, increasing the resolution leads to larger gaseous 
disks, but in some cases the total mass in gas increases while 
it decreases in others. These large variations are at least 
partly due to the fact that some galaxies (e.g., R-AGN) have 
almost no gas left at z = 0, and therefore even small changes 
can lead to large fractional variations. The properties of the 
gaseous component seem the most vulnerable to numerical 
resolution effects and therefore caution must be exercised in 
their interpretation. 



4 SUMMARY AND CONCLUSIONS 

The Aquila Project compares the results of simulations of 
galaxy formation within a Milky Way-sized ACDM halo. 
We use 9 gas-dynamical codes developed and run indepen- 
dently by different groups adopting in each case their pre- 
ferred implementation of radiative cooling, star formation, 
and feedback. The codes include SPH-based, grid-based, 
and moving-mesh techniques; all include supernova feedback 
in thermal and /or kinetic form and primordial or metal- 
dependent cooling functions. Two of the codes (gadget and 
RAMSES) were run three times with three different sub-grid 
physics for a total of 13 different simulations. In addition, 
each code was run at two different resolution levels in or- 
der to investigate numerical convergence. All runs share the 
same initial conditions and are analyzed with a set of con- 
sistent analysis tools. Our main conclusions may be summa- 
rized as follows. 

Galaxy Morphology. At « = 0, simulated galaxies ex- 
hibit complex morphologies, with spheroids, disks, and bars 
of varying importance. Morphology shows no obvious depen- 
dence on the hydrodynamical method adopted or on numer- 
ical resolution, and seems to be mostly related to how and 
when gas is accreted and transformed into stars. The best in- 
dicator of the presence of a disk seems to be the median age 
of the stars; the later stars form the more prevalent the disk 
component is. This suggests that, to be successful at forming 
disks, codes must be able to preempt the early transforma- 
tion of gas into stars while at the same time promoting the 
accretion and retention of late-accreting, high-angular mo- 
mentum gas. 

Stellar Mass. Despite the common halo-assembly his- 
tory, we find large code-to-code variations in the final mass of 



© 0000 RAS, MNRAS 000, 000-000 



18 Scannapieco, Wadepuhl, Parry, Navarro et al. 



simulated galaxies. The stellar mass ranges from4xlO^°M0 
to ~ 3 X 10^^ Mq, depending largely on the adopted strength 
of the feedback. There is also an indication that the nu- 
merical method may play a role: arepo is able to form al- 
most twice as many stars as G3 although they both adopt 
the same sub-grid physics. Most simulated galaxies are too 
massive compared with theoretical expectations based on 
abundance-matching considerations. The median stellar age 
also correlates with galaxy mass, indicating that models that 
favor late star formation (as needed to form a disk) do so at 
the expense of allowing too many stars to form overall. 

Rotation Curves . All simulated galaxies have rotation 
speeds in excess of what is expected from the Tully-Fisher 
relation of late-type spirals. The disagreement worsens as 
the stellar mass of the simulated galaxy increases, both 
for galaxies with and without a well-defined stellar disk. 
At the high mass end, simulated galaxies have unrealisti- 
cally sharply-peaked, strongly declining rotation curves. Al- 
though reasonably fiat rotation curves are obtained at low 
Mjteiiar, the TuUy-Fisher zero-point offset persists for these 
systems. 

Galaxy Size. The difficulties matching the Tully- 
Fisher relation are due to the fact that most simulated galax- 
ies have stellar half-mass radii much smaller than expected 
from observation given their stellar mass. This is especially 
true for galaxies where inefficient feedback allows the for- 
mation of an overly massive galaxy; these systems are, quite 
unrealistically, smaller even than dense early-type galaxies. 
Galaxies where feedback is able to limit the stellar mass 
to more acceptable levels are also more concentrated than 
typical spirals, highlighting the difficulty that all codes face 
to prevent too many low-angular momentum baryons from 
assembling into the galaxy. 

Star Formation History. The excess of low-angular 
momentum baryons alluded to above may be traced to the 
inability of feedback schemes to prevent large bursts of star 
formation at early times. This is clearly seen in the stel- 
lar age distribution, which shows that star formation peaks 
at a ~ 4 and declines steadily afterwards. The same trend 
holds for essentially all runs, albeit modulated by the partic- 
ular implementation of feedback adopted by each individual 
code. Indeed, essentially all models allow more than half of 
the stars to form in the first ^ 3 Gyrs of evolution. The 
relative insignificance of star formation in recent times com- 
pared to the intensity of the early burst seems to be at the 
root of many of the difficulties that simulated galaxies have 
in matching observation. 

Gas Component . The properties of the gaseous compo- 
nent of the galaxy at z = show even wider code-to-code 
variations than the stars, since it is continuously resupplied 
by accretion, depleted by star formation, and dislodged by 
feedback- driven winds. This results in large short-lived fluc- 
tuations that lead to poor numerical convergence and large 
code-to-code scatter. Most simulated galaxies have gaseous 
disks with sizes that compare favorably with observation, 
although their gas fractions are systematically lower than 
those of star-forming spirals of comparable mass. 

Numerical Convergence . We have investigated numer- 
ical convergence by comparing the results at two different 
numerical resolutions, spanning a range of ~ 2 and ~ 8 in 
spatial and mass resolution, respectively. Although numeri- 
cal convergence is not particularly good for any of the codes, 



reasonably-good convergence is found for the properties of 
the stellar component, such as total mass and median age. 
Less well-converged are the internal properties of the galaxy, 
such as the half-mass radius, or the fraction of stars in a 
rotationally-supported disk. For the same reasons cited in 
the previous item, the properties of the gas are the ones 
that are most poorly reproduced at the two different reso- 
lutions. There is also indication that feedback schemes that 
are more effective at limiting the stellar mass of the galaxy 
(such as through the inclusion of AGN-related feedback in 
addition to supernovae) converge less well than other imple- 
mentations. Overall, the variations introduced by resolution 
are small compared to code-to-code variations, which leads 
us to conclude that none of the trends highlighted above are 
driven primarily by resolution. 

SPH vs AMR vs Moving-Mesh. Our results suggest that 
numerical hydrodynamics techniques have some influence on 
the outcome of a simulation. This is most clearly demon- 
strated by comparing the results of G3 and AREPO which, 
despite adopting the same sub-grid physics modules, yield 
galaxies that differ by almost a factor of 2 in stellar mass. 
The AMR technique also yields large stellar masses (when 
similarly inefficient feedback is adopted, as in run R), lend- 
ing support to the view that standard SPH-based codes may 
underestimate the total amount of gas that can cool and be- 
come available for star formation at least in the weak feed- 
back regime. On the other hand, the galaxies formed by R, 
G3, and arepo are unrealistically massive and concentrated, 
so large modifications to the feedback implementation of 
these codes are needed in order to bring them into agree- 
ment with observation. These changes may overwhelm the 
method-induced differences; for example, R makes a promi- 
nent disk while R-AGN has 5 times fewer stars and no disk. 
It is therefore unclear at this point whether the shortcom- 
ings of SPH are actually significant compared with the un- 
certainties involved in designing a star formation/feedback 
scheme that can yield realistic galaxies, although it is ob- 
viously desirable to avoid numerical inaccuracies as far as 
possible. Further investigation of this question is needed to 
clarify this issue. 

Aside from these considerations, perhaps the main re- 
sult of the Aquila Project is that, despite the large spread in 
properties spanned by the simulated galaxies, none of them 
has properties fully consistent with theoretical expectations 
or observational constraints in terms of mass, size, gas con- 
tent, and morphology. Despite this apparent failure, we be- 
lieve that the Aquila Project nevertheless yields interesting 
clues to guide how codes might be modified to yield realis- 
tic galaxies. For example, the need (i) to control effectively 
the overcooling of baryons; (ii) to curtail the early burst 
in star forming activity; and (iii) to promote the accretion 
and retention of the late-accreting, high-angular momentum 
baryons needed to form disks similar to those of normal spi- 
rals are of general applicability to all codes. 

There seems to bo little predictive power at this point in 
state-of-the-art simulations of galaxy formation; these seem 
best suited to the identification of the role and importance 
of various mechanisms rather than to the detailed modeling 
of individual systems. It may be argued that the strength 
of this conclusion depends on whether the parent halo of 
the Aquila runs (Aq-C) is truly destined to harbor a disk 
galaxy and that there is no hard proof for this. Further, the 
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possibility that Aq-C might be an unrepresentative outUer 
should also be considered, as suggested by the l-GALAXIES 
semi-analytic model (see, e.g., Fig. O. 

These objections may be addressed when simulations 
are able to follow a statistically-significant number of galax- 
ies in a cosmologically-representative volume. We might not 
know what kind of galaxy inhabits an individual halo, but we 
do know what the population of galaxies looks like. Evolving 
a region large enough to contain at least a few dozen Milky 
Way-sized galaxies at the resolution achieved here seems like 
a natural next step, and one that should be achievable in the 
not too-distant future. 

Finally, the complexity of the problem suggests that 
the best approach to improving galaxy formation simula- 
tions may be one where multiple numerical alternatives are 
developed and explored simultaneously and independently, 
provided that they are periodically contrasted in controlled 
experiments such as the one presented here. Given the in- 
tricacy of the task and in the absence of a clear front- 
runner, the diversity of numerical techniques presently avail- 
able might actually turn out to be a strength rather than a 
shortcoming. 
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Table 3. Properties of simul ated galaxies at z = 0, for the leve l 5 (upper row) and level 6 (lower row) resolution simulations. We also list the prediction of the semi-analytic models 
GALFORM jCooper et alj[20ich and L-GALAXIES jGuo et alj|201 j) . and of the dark matter-only simulation Aq-C-4 of the Aquarius Project. 
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0.10 
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0.35 
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APPENDIX A: DESCRIPTION OF THE CODES 

In this Appendix we summarize the main properties of the 
various simulation codes and implemented physics. These 
succinct descriptions have been provided by each of the in- 
dividual groups participating in the Aquila Project. As ex- 
plained below, the different models use a variety of imple- 
mentations of the processes of star formation and feedback, 
resulting in the star formation rate surface densities shown 
in Fig. 

Gadget-3 models (G3, G3-BH, G3-CR) 

The simulations G3, G3-BH and G3-CR are based on 
GADGE TS. This is a s ignificantly updated version of GAD- 
GET2 (|Springel|[200H ). a fully cosmological code based on 
smoothed particle hydrodynamics (SPH) and on the Tree- 
PM method to evaluate gravitational forces. All three of 
th ese simulations use th e star formation model introduced 
bv lSpringel fc Hernguistl ll2003l) . G3-BH includes AGN feed- 



back following 
model includes 



Springel et al 



(|2005al '). while the G3-CR 

addition, energy deposition by cosmic 

rays (see Jubelgas et al. 20081. and re ferences therein). 

The lSpringel fc Hernguistl (|2003h star formation model 
uses a primordial cooling function following iKatz et al.l 
l|l996l) and includ es a uniform UV backgro und, based on an 
updated version of lHaardt fc Madaul l| 19961) . that is switched 
on at 2 = 6. Stars are formed stochastically in dense regions 
assuming a Salpeter IMF. This model pictures each gas el- 
ement as a two-phase mixture of hot and cold gas in ap- 
proximate pressure equilibrium. Cold gas is converted into 
"star" particles on a characteristic timescale t*. Assuming 
that a fraction /3 of the stellar population represented by 
each star particle is so short-lived that they explode as su- 
pernovae instantly, their (metal-enriched) mass is fed back 
to the interstellar medium. The star formation rate is given 

by 



dp, 
dt 



(Al) 



where pc is the mean mass density of gas in the cold phase, 
and is a density-dependent star formation timescale. The 
energy feedback of supernovae explosions is used to heat the 
ambient hot gas as well at to evaporate cold clouds, leading 
to a self-regulating cycle for star-forming gas. The net effect 
is that the dynamics of the multi-pha se medium is governed 
by an effective equation-of-state (see ISpringel fc Hernguistl 
I2OO3I ). 

The simulations G3-BH and G3-CR both include ther- 
mal AGN feedback associated with a Bondi-Hoyle-Lyttelton 
parameterization of the (spatially unresolved) gas accretion 
onto a central supermassive black hole (see iSpringcl et al.;l 
l2005al . for a detailed description) . The thermal feedback en- 
ergy is parameterized by Eiccd = e/erMsHC^, where is 
the radiative efficiency of the black hole and e/ encodes the 
feedback efficiency. 

In G3-CR the distribution function of cosmic rays (CR) 
is modeled as a power law. 



d^N 
dp dV 



(A2) 



in momentum space. The non-thermal pressure of this cos- 
mic ray population is then given by 



CmpC2 
>x;r = o 1 

6 1+92 



2 3 



(A3) 



which is added to the ordinary gas pressure jjubelgas et al.l 
2008). Here B„{a,b) denotes incomplete Beta functions and 
Q is the assumed power-law slope of the CR particles. The 
energy injected by supernova explosions into the CR pop- 
ulation is given as AesN = CsNesNJTi^Af where ("sn is the 
is the fraction of the supernova energy that first appears in 
cosmic rays. 

The numerical parameters used for the Aquila Project 
simulations are as follows: a star formation timescale oft* — 
2.1 Gyr and a supernova energy of Esn = 10""^ erg. The G3- 
BH model adopts = 0.1 and e/ = 0.05. The G3-CR run 
uses a cosmic ray generation efficiency of (^sn = 0.3 and a 
spectral index a = 2.5 

The CS model (G3-CS) 

The CS model is a GADGET3-based code that in- 
cludes stochastic star formation, chemical enrichment and 
supernova (thermal) feedback from Type II and Type la 
SN explosions, a multi-phase model for the gas compo- 
nent and metal-dependent cooling. Full detail s of t he im- 
plementation are given in IScannapieco et all (2005) and 
IScannapieco et al] l|2006l ) and previous applications to cos- 
m ological galaxy formation and d i sk formation can be found 
in IScannapieco et al.l (|2008l . l2009l . l20ld . |2011^ . 

The model includes a UV background field which 
tur ns on at ^: = 6 a nd follows the formulation 
of iHaardt fc Madaul d 199(3) and metal - depen dent cool- 
ing functions from Sutherland fc Dopital (|l993l ) are used. 
Star formation occurs stochastically in gas particles 
with high density (n > nth ) and in a convergent 
fiow (|Springel fc Hernguistl liooj ). following the Schmidt- 
Kennicutt law: 



dp, 
~dt 



^dyn 



(A4) 



where p» and pgas are the stellar and gas densities, c« is a star 
formation efficiency and tdyn is the dynamical time of the 
gas. Stars return metals and energy during supernova Type 
II and Type la explosions. Both types of SNe eject the same 
amount of energy _Esn to the interstellar medium, but they 
have different chemical yields, explosion times and rates. SN 
energy is dumped in given fractions to the cold (ec) and hot 
(1 — Ec) neighbors of exploding stars. Energy deposition into 
hot neighbors occurs at the time of explosion; however, for 
cold neighbors energy from successive explosions is instead 
accumulated and deposited only after a time-delay which de- 
pends on the local conditions of the cold and hot gas phases. 
In this way, artificial loss of SN energy in high-density re- 
gions is prevented. Star particles are treated as single stellar 
populations with a Salpeter initial mass function. 

We introduce a multiphase gas model which allows gas 
in both dense and diffuse phases to coexist in the same spa- 
tial region. In our model, SPH particles of a given physi- 
cal state (density, entropy) ignore neighboring particles that 
have a much lower (a factor of 50) entropy. This scheme also 
makes the deposition of SN energy more efficient. 

The input parameters used for the simulations ana- 
lyzed in this paper are: nth = 0.03 cm~^, c, — 0.1, 
i?SN ~ 0.7 X 10^^ erg and ec ~ 0.5. Finally, the input pa- 
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Figure Al. Cold (T < 10^ K) gas surface density versus star formation rate per unit area in tlie simulated galaxies, measured face-on 
and averaged within a cylinder of radius equal to the (projected) half- mass radius of the cold gas (see the right panel of Fig. [9J. Each 
set of symbols indicate results for level-5 simulations aX, z = 2,1.5,1,0.5 and 0, with the size of the symbols increasing with decreasing 
redshift. The das hed line indicates, for reference, the Schmidt-Kennicutt law for nearby "normal" and "star-bursting" disks, as given by 
iKennicuttI l|l99Sl 'l. 



ra meters for the chemical m odel are identical to those used 
in lScannapieco et al.l (|2009l '). 

The TO model (G3-T0) 

The TO model is described in lOkamoto et all (|2010l ') 
and is based on an early version of the GADGETS code. 
It incorporates metal-dependent cooling, star formation, 
thermal and kinetic feedback from supernovae (SNe) and 
enrichment by AGB stars and SNe. Examples of its ap- 
plicat ion in cosmologica l galaj c y formation sim u lation s in- 
cludelOkamoto fc Frenl3 (|2009l ). lOkamoto et~al] l|201(]| ) and 
iParrv et al.l (|201ll ). 

Photo -heating and coohng are implemented as de- 
scribed in |w lersma et al.l (|2009al ). including contributions 
from eleven elements in the presence of a spatially uni- 
form, ti me evolving UV background of the form calcu- 
lated bv lHaardt fc Madaul ([2001). Gas above a density nth 
forms stars at a rate normalized to reproduce the Schmidt- 
Kennicutt law. Star part i cles re present simple stellar pop- 
ulatio ns with a | Chabrieil (|2003| ) initial mass function. Fol- 
lowing |w iersma et al. I (|2009bl ). energy, mass and metals are 
returned to the ISM by AGB stars, type la and type II 
SNe on timescales appropriate for the age and metallicity of 
the s tellar population, with yiel ds and stellar li fetimes taken 
from lPortinari et all (|l998l ) and lMarigol (|200ll ). 

The unresolved interstellar medium is modeled using 
a sub-grid prescription. Each gas particle with n > nth is 
treated as a series of cold clouds, with an empirically moti- 
vated mass function, embedded in a hot ambient medium. 
The two phases exchange mass through thermal instability 



and cloud evaporation. Each cloud has a star formation rate 
that is inversely proportional to its dynamical time. The 
hot phase is supported by imposing a minimum pressure of 
Pmin oc p^'^. Type la SNe increase the thermal energy of the 
gas around star forming regions, whilst type II SNe are as- 
sumed to drive large-scale winds. The wind speed is chosen 
to be proportional to the halo circular velocity, using the lo- 
cal dark matter velocity dispersion as a proxy (v — oigdm)- 
The mass loading then follows by requiring conservation of 
all of the available SNe energy. Wind particles are decoupled 
from hydrodynamic forces for a short time, to allow them to 
escape the star forming region and ensure that the specified 
mass loading and wind velocity are not modified by viscous 
drag from the ISM. 

The density required for star formation is set at 
nth =0.1 cm~^ in the level 6 simulation and 0.4 cm~"^ in 
the level 5 simulation. For the wind speed parameter, q = 5 
is used, which has been shown to produce a good match to 
the Milky Way satellite luminosity function (|Okamoto et al.l 
I2OIOI : iParrv et al.ll201ll ) . 

The GIMIC model (G3-GIMIC) 

The GIMIC model is a GADGExS-based code that in- 
cludes metal dependent cooling on an element-by-element 
basis in the presence of a UV background, star formation 
and supernovae (SNe) driven winds, as well as mass and 
metal recycling by AGB stars, type la and type II SNe. 
GIMIC is id entical to model M ILL of the OWLS suite of 
simulations jSchave et al. 20101). A full descr iption of the 
model can be found in ICrain et all (|2009l ) and lSchave et al.l 
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(l201(t) and further applications of thi s model can be found 
in lCrain et al] l)2010l).ICui et all ll201lll,lDeason et al.l (|201ll ') 



and lFont et alj (|201ir )T lMcCarthv et all (|201ll ). 

The model includes a spatially uniform, time 
evolving UV ba c kgrou nd following the formulation of 
iHaardt fc Madaul (|200ll 'l. Hydrogen reionises at z; = 9, 
Helium H at 2; = 3.5. Radiative cooling and heat- 
ing processes are implemented on an element-by-element 
basis, using in t erpol ation tables compute d with cloudy 



Ferland et al.l 1 19981 ). as described by IWiersma et al.l 
2009al ). Th e star formation rate (S F R) p rescription, de- 



scribed by ISchave fc Dalla Vecchial (|2008l ). IS pressure- 
dependent and enforces a local Schmidt-Kennicutt law and 
requires only that the slope and normalization of the ob- 
served relation are specified as input parameters. Gas par- 
ticles are converted to star particles stochastically, with a 
probability that depends on their associated SFR. Each star 
particle r e prese nts a simple stellar population (SSP) with a 
IChabrieJ (|2003l ') initial mass function and inherits the ele- 
mental abundances of its parent gas particle. The chemo- 
dynamical evolution of SSPs and the associated recycling 
of heavy elements into surrounding gas, is followed on an 
element-by-element basis and includes contributions from 
AGB stars a nd both Type l a and Type II supernovae, as 
described bv lWiersma et al.1 (12009^ 1. 

GIMIC appeals to a phenomenological treatment to 
model the energetic feedback resulting from stellar evolution 
and supernovae. As the simulation lacks both the physics 
and the resolution to model the multiphase ISM, an effec- 
tive equation of state is imposed onto gas particles that are 
sufficiently dense (nn > 0.1 cm"^) an d cold (T < 1 0^ K) to 
be subject to gravitational instability (|Schavell2004l '). The ef- 
fective equation of state, P — Kp^l'^ , is chosen to ensure that 
the Jeans mass and the ratio between the Jeans length and 
the SPH smoothing length are independent of the gas den- 
sity, thus preventing spurious fragmentation due to a lack of 
numerical resolution (|Schave fc Dalla Vecchiall2008l ') . SN en- 
ergy is assumed to drive gala:xy-wide outflows, therefore gas 
particles neighboring star-forming regions are given a ran- 
domly orientated velocity kick of 600 kms~^. A probabilistic 
scheme ensures that, on average, the mass put into the wind 
is a factor of four times the amount of stars formed. Unlike 
many similar schemes, these particles are not temporarily 
dec oupled from hydrodynamic for ces (for further discussion, 
see iDalla Vecchia fc Sch ave"2008^) . 

Values for a ll model parameters can be found in 
ICrain et al.1 (|2009l ). but to summarize: the initial wind mass 
loading and velocity are set to 4 and 600 kms^^ respectively 
(with 10^^ ergs per SN), the star formation threshold to 
n > 0.1 cm~^ and the input star formation law to have slope 
1.4 and normalization coefficient ~ 1.5 x lO~'*M0yr~^kpc~^. 

The MM model (G3-MM) 

The MM model (or M Ulti Phase Part i cle In tegrator, 
MUPPI), fully described in iMurante et al.1 l|2010l '). is also 
implemented within GADGETS. We include radiative cooling 
of a gas with primordial composition; he ating from a uni- 
form U V background of the form given bv lHaardt fc MadarJ 
l|l996ll . turned on at a redshift 2 = 6; star formation and 
stellar feedback with the multi-phase model described be- 
low. We assume a Salpeter IMF and the Instantaneous Re- 
cycling Approximation to treat gas restoration and SN en- 



ergy feedback. No chemical enrichment is included in the 
present version of the code. No treatment of accreting black 
holes or cosmic ray feedback is implemented. 

In the MUPPI sub-resolution model of star formation 
and stellar feedback, gas particles at relatively high density 
(n > nth) and low temperature (T < Tth) are treated as 
a multi-phase system, made up by cold and hot gas phases 
coexisting in pressure equilibrium, and a stellar component. 
Hot phase has low mass fraction but high filling factor, so its 
cooling time is much longer than that computed for the aver- 
age particle density, and thus thermal energy is not quickly 
dissipated as soon as it is injected. Regarding the cold phase, 
a part of it is assumed to be in mol ecular form. Blitz fc 
Rosolowsky (|Blitz fc Rosolowskvl[2006h found a phenomeno- 
logical correlation between the ratio of molecular and HI gas 
surface densities and the external disk pressure. Inspired by 
this result, we compute the fraction of cold gas in molecular 
form as a function of hydrodynamical pressure P: 



/mol — 



1 + 



(A5) 



Mass flows among the components as follows. Cooling 
of the hot phase feeds the cold gas phase. Stars form from 
the molecular cold gas: 



tdy 



(A6) 



where fdyn is the dynamical ti me of the cold phase at the 
onset of a multi-phase cycle (see lMurante et al1l201Gl for de- 
tails) and /* is the fraction of molecular cloud transformed 
into stars per dynamical time. Our prescription for star for- 
mation is not based on imposing a Schmidt-Kennicutt re- 
lation, which is instead naturally produced by our model 
(see Monaco et al. 2012). A fraction /ov of the star-forming 
gas is evaporated back to the hot phase, while a fraction 
/rc of stellar mass is instantaneously restored to the hot 
phase. The code also tracks the hot phase thermal energy. 
This gained or lost through hydrodynamics, lost by cool- 
ing and gained from type II SN explosions: a small frac- 
tion /in of their energy is deposited into the hot phase of 
the same star-forming gas particle, while a more significant 
fraction /out is distributed to neighboring particles, prefer- 
entially along the "least resistance path" defined as a cone 
of semi-aperture 6 and anti-aligned with the density gradi- 
ent. Finally, to mimic the destruction of molecular clouds a 
multi- phase particle returns single- phase after a time tciock, 
scaled with the dynamical time. 

A system of ordinary differential equations evolves the 
mass and energy flows described above. This is solved on- 
the-fiy within each SPH time-step with a Runge-Kutta in- 
tegrator with adaptive time-step. This dynamical system 
has an intrinsically runaway behavior: energy from SNe in- 
creases gas pressure, which in turn increases SFR through 
the higher molecular fraction. This runaway is stabilized by 
the hydrodynamic response of gas: when a particle receives 
enough energy, it expands thereby decreasing its pressure. 
The equations MM model solves are similar to those used in 
the star formation model by Springel & Hernquist (2003); 
the main differences are that in MM model no equilibrium 
solution is assumed, and the effect of hydrodynamics on the 
multi-phase gas is explicitly taken into account. 

MUPPI produces reservoirs of "virtual" stars, that are 
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transformed into star particles using the stochastic algo- 
rithm of Springel & Hernquist (2003). Each gas particle pro- 
duce up to 4 generations of star particles. 

We use d here the sam e "sta ndard" set of parameters 
described in lMurante et all (|2010t ). namely: (i) star forma- 
tion efficiency /* = 0.02, amount of molecular gas which 
is converted into stars; (ii) Po ~ 35000 K cm~^, pres- 
sure normalization for the Blitz & Rosolowsky relation; (m) 
Tc = lOOOK, temperature of the cold phase; (iv) font = 0.3, 
fraction of SNe energy given to neighboring gas particles; 
(v) fin = 0.02, fraction given to the hot gas of the par- 
ticle itself; (vi) f^v = 0.1, fraction of cold gas mass evap- 
orated by SNe; (vii) Esn = 10^""^ erg, SN energy; (vii) 
6 = 60, semi-aperture of the cone determining the neigh- 
boring gas particles which receive energy; (viii) nth = 0.01 
cm~^, density threshold for entering the multi-phase stage; 
(ix) Tth = 50000K, temperature threshold for entering the 
multi-phase stage; tciock = 2fdyn, time after which a particle 
exits the multi-phase stage. 

The CK model (G3-CK) 

The CK model is a GADGET3-based code that includes 
star formation, chemical enrichment, and (thermal) feedback 
from stellar winds, core-collapse supernovae. Type la su- 
pernovae, and asympto tic giant b r anch ( AGB) stars. The 
details are described inlKobavashil ||2004|) . [Kobavashi et al.l 
(|2007tl , and lKobavashi fc Nakasatol l|201 j ). 

The UV backgrou nd radiation is included with 
iHaardt fc Madaul (| 19961 ) from z — &. Radiative cooling 
is computed with the metal-dependent cooling functions 
l|Sutherland fc Dopital fToQa^ as a function of [Fe/H]. [O/Fe] 
is fixed with the observed [O/Fe]- [Fe/H] relation in the solar 
neighborhood. The star formation criteria are 1) converging 
flow, 2) rapid cooling, and 3) Jeans unstable. The star forma- 
tion rate (SFR) is determined from the Schmidt-Kennicutt 
law (Eq. IA4|) . If a gas particle satisfies the star formation 
criteria, a part of the mass of the gas particle turns into a 
new star particle. We then treat the star particle as a simple 
stellar population and calculate the evolution of the stellar 
population every time-step. The masses of the stars asso- 
ciated with each star particle are distributed according to 
an initial mass function. We adopt the Salpeter IMF that is 
invariant to time and metallicity with a slope x = 1.35 for 
0.1 — 120 M(;), to be consistent w ith the Galactic chemical 
evolution (Kobavashi et al1l2006l ). 

For the feedback of energy and heavy elements, we do 
not adopt the instantaneous recycling approximation. In- 
stead, via stellar winds, core-collapse supernovae. Type la 
supernovae, and AGB stars, thermal energy and heavy el- 
ements are ejected from an evolved star particle as a func- 
tion of time, and distributed to a constant number A^fb of 
surrounding gas particles. Among core-collapse supernovae, 
we include the effect of hypernovae, which are observation- 
ally known to produce more than ten times larger explo- 
sion energy and iron than normal Type II supernovae. We 
adopt the metal-dependent efficiency of hypernovae (ehn = 
0.5,0.5,0.4,0.01, and 0.01 for Z = 0,0.001,0.004,0.02, and 
0.05) for the initial masses of M > 2OM0, to be consis- 
tent wit h the observed [Zn/Fe ] rati os in the Milky- Way 
Galaxy (|Kobavashi fc Nakasatd I2OIII ). The ejected energy 
(1 — 40 X 10^^ erg) and nu cleosynthesis yields are taken 
from iKobavashi et al.l (|2006l ) as a function of progenitor 



mass and metallicity. With hypernovae, cosmological sim- 
ulatio ns give a better agre ement with observed cosmic 
SFRs (|Kobavashi et al.|[2007l '). For Type la Supernovae, we 
use our single-d egenerat e model with the meta llicity effect 
IKobavashi et al.. .199S ; IKobavashi fc Nomotd ^009). The 
lifetimes span a range of 0.1 — 20 Gyr as a function of pro- 
genitor metallicity, which is consistent with the observed su- 
pernova rates and with the observed [a/Fe] relations in the 
Milky- Way Galaxy. The ejected energy is 1.3 x 10^^ erg per 
explosion. ^From stellar winds, ~ 0.2 x 10^^ erg is ejected 
depending on metallicity for the stars with > 8Mq. The 
adopted nucleosynthesis yiel ds of supernov ae and AGB stars 
are summarized in Kobavas hi et al.l (|201ll ) . 

The input parameters used for the simulations analyzed 
in this paper are: the star formation timescale of c, = 0.02 
and the feedback number of A^fb = 64. 

The Gasoline model ( GAS) 

The GAS model uses the SPH code g asoline, which is 
described in detail in IWadslev et al] (|2004l '). It includes a UV 
background heating, metal cooling, star formation, thermal 
stellar feedback and chemical enrichment from SNII, SNIa, 
and mass return from stellar winds. 

GASOLINE contains metal cooling based on radiativ e 
transfer found in CLOUDY as described in IShen et al.ll|2010l ). 
The CLOUDY cooling was calculated using an external UV 
radiation field starting at z = 8.9. Star formation is cal- 
culated and supernova feedback is i mplemented us i ng th e 
blastwave formalism as described in IStinson et al.l (|2006h . 
Recent examples of simulations that use similar physics 
includ e iGov crnato ct al. (2007, 2009, 2010), Stinso n et all 
(|2010l ). lBrooks et al] (|201ll ) and lGuedes et al.l (|201li ). 

Star formatio n is calculat ed using the commonly used 
recipe described in lKatzll|l992h . Stars form from gas below a 
maximum temperature of 15,000 K and above a density of 1 
cm^'^ with an effi ciency of 5%. However, the h igh resolution 
runs p resented in lGovernato et al] (|2010l ) and lGuedes et al] 
(I2OIII ) used a higher threshold for star formation (100 and 
5 cm~^ respectively) that leads to more efficient gas out- 
flows than the SF recipe adopted for the Aquila simula- 
tion. The star particles are treate d as single ste l lar po pu- 
lations using the IMF described in iKroupa et all (|l993l '). In 
this context, ejecta from both Type II and type la super- 
novae are considered. These supernovae feed both energy 
and metals back into the interstellar medium gas surround- 
ing the region where they formed. Type II supernovae de- 
posit 7 x 10^" erg of energy into the surrounding interstellar 
medium. Since this gas is dense, it would be quickly radi- 
ated away due to the efficient cooling. For this reason, cool- 
ing is disabled for particles inside the blast region rcso 
10 
10 
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6.85 ciO. 32^0.34 i 

iisi no 1 



04 
-0.70 



pc for t he length of time tcso = 
yr given in lMcKee fc Ostrikeij (|l977t) . 
In summary, the input parameters used for the simula- 



tions analyzed in this paper are: nth 
and Esn = 0.7 x 10^^ erg. 



1.0 cm-^ c. = 0.05, 



The Ramses models (R, R-LSFE, R-AGN) 

The RAMSES code is an Eulerian Adaptive Mesh Refine- 
ment code that uses the Particle Mesh techniques for the N 
body part (stars and dark matter) and a shock-capturing, 
unsplit second-order MUSCL scheme for the fluid compo- 
nent. For the latter, we use the MinMod slope limiter and the 
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HLLC Riemann solver, ens uring both stab ility and proper 
capturing of discontinuities (|Tevssieill2002l ). 

Star formation is implemented stochastically using a 
Schmidt law, with density threshold for star formation held 
fixed at nth ~ 0.1 cm~^ and an efficiency parameter chosen 
between 1% (R-LSFE) and 5% (R, R-AGN). Stellar feedback 
is modeled using a thermal dump of 10^^ erg per supernova, 
assuming a Salpeter IMF. Cooling is performed using a tab- 
ulated cooling function depending on gas metallicity, the lat- 
ter being modeled self-consistently as a additional scalar hy- 
dro variable an d initiated during supernovae explosions with 
a yie ld y=10% (|Rasera fc Tevssieij|2006l : iDubois fc Tevssieij 
l2008l ). 

The grid is refined following a quasi-Lagrangian strat- 
egy, each cell being refined if the number of dark matter 
particles exceed 8, or if the baryonic mass (star + gas) ex- 
ceeds 8 times the initial mass resolution set by the Aquila 
initial conditions. In order to avoid catastrophic refinement 
at early times, we carefully trigger new levels so that the 
maximum level levelmax is opened only at expansion fac- 
tor a — 0.8, the previous one at a = 0.4, levelmax— 2 at 
a = 0.2 etc. This ensures that the resolution of the grid in 
physical units remains roughly constant (although in a dis- 
crete, stepwise sense). For the level 6 simulation, we have 
set levelmax=17 (cell size is 1 kpc physical) and for the 
level 5 simulation, levelmax=19 (cell size is 261 pc physi- 
cal). This corresponds also to the maximum levels reached 
by a pure dark matter simulation with the same number of 
dark matter particles for each case, minimizing spurious ef- 
fects due to two-body relaxation. For the R-AGN simulation 
o nly, AGN feedback has b een i mplemented us i ng th e model 
of iBooth fc Schavd (|2010l ): see lTevssier et al.l (|201ll ') for de- 
tails. A sphere of size 4 cells is defined around each sink 
particle defining the SMBH. This spherical region is used 
to determine the accretion rate on the SMBH, but also to 
spread on the grid the AGN feedback energy, using only a 
pure thermal dump. 

To summarize, the RAMSES simulations use a star for- 
mation density threshold of nth = 0.1 cm~^, an energy per 
SN of lO^'^ ergs and a star formation efficiency of 5% (R and 
R-AGN) or 1% (R-LSFE). 

The Arepo model (AREPO) 

The AREPO code (ISpringellliblOal ) is a novel pseudo- 
Lagrangian hydrodynamical code that works with an un- 
structured, fully dynamic and adaptive mesh. The m esh is 
defined as the Voronoi tessellation (e.g. IOkabe|[200Cl ) of the 
simulation volume generated by a set of mesh-generating 
points. The hydrodynamics is calculated with a second-order 
accurate finite volume approach on this mesh, based on the 
MUSCL-Hancock scheme that is also widely employed in 
ordinary Eulerian mesh codes. This involves a spatial recon- 
struction step and flux estimates at all cell faces by solving 
Riemann problems at the interfaces. The very important 
new ingredient in the scheme is that the mesh-generating 
points are allowed to move freely during the calculation. In 
particular, they can be moved with the local fluid veloc- 
ity such that the mesh dynamically follows the fluid motion 
without showing any problematic mesh twisting effects. In 
this mode, AREPO minimizes advection errors in the hydro- 
dynamics and produces Galilean-invariant results that help 
to avoid accuracy problems with high velocity cold flows that 
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Figure Bl. Evolution of the physical gravitational softening in 
the different models. 

can occur in ordinary mesh codes. A more detailed descrip- 
tion and an investigation of numerous test problems that 
demonstrate the high accuracy of the method can be found 
in Springol (2010a). 

Even though the hydrodynamics is solved completely 
differently in AREPO than in GADGETS, the Lagrangian 
character of both codes makes it possible to implement the 
physics of star formation and feedback in very similar ways 
in both codes. In fact, the AREPO simulations analyzed here 
implement exactly the same sub-grid model for dense gas as 
well as the same supernova feedback recipe as our default 
GADGET run G3. Since also the gravitational solver for both 
codes is identical, any differences found in the results hence 
reffect changes due to the numerical treatment of hydrody- 
namics alone. 



APPENDIX B: GRAVITATIONAL SOFTENING 

Our 13 simulations have used a variety of choices for the 
evolution of the gravitational softening. This evolution is 
governed by two parameters: ^ax and el"". The former di- 
vides the period where the softening changes from being 
fixed in comoving coordinates [z > Zflx) to being fixed in 
physical coordinates {z < zax). The different simulations 
have various choices for Zhx, as shown in shown in Tabled 
The second parameter, e^^", is the value of the gravitational 
softening at the present time, and also varies slightly for our 
13 simulation (Table [2]). Fig. IB II shows the evolution of the 
gravitational softening in physical coordinates for our simu- 
lations. 



APPENDIX C: DISK ORIENTATION 

Fig. [Cl] illustrates, for our level 5 and level 6 runs, the spin 
of the stellar component of the simulated galaxies at z — 0. 
The two panels show, in a 3D rendering, the specific angu- 
lar momentum vector of all stars in each galaxy, normal- 
ized to the maximum among all simulations. Clearly, disks 
are not all aligned in the same direction, nor is the spe- 
cific angular momentum the same. This is not surprising, 
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given the wide range of stellar masses spanned by the differ- 
ent simulations. In spite of this, the orientation is actually 
similar for some simulations. As explained in Section 12.41 
for orientation-dependent diagnostics we have rotated each 
simulated galaxy to a new coordinate system where the an- 
gular momentum vector of its stellar component coincides 
with the z direction. 
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Figure CI. Disk orientation of the level 5 (left-hand panel) and level 6 (right-hand panel) resolution simulations. The vectors indicate 
the orientation of the angular momentum of galactic stars, and are normalized to the maximum among all simulations. 
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The MNRAS K-TJ^X 2e class file is implemented by placing the command 
\documentclass [optionsl ■[nin2e}- 
at the start of the document. 

The various option commands listed in the Monthly Notices MJ^X style guide for authors can still be used. The following 
additional options exist. 
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fonts (eurmm); characters ^ and appear slanted only on systems that have the AMS series A fonts (msamxx). On systems 
that do not have these fonts, the standard forms of the characters appear in the printout; however, they should be correct in 
the final typeset paper if the correct I^T]gK commands have been used. 
Q ■ (ii) usedcolumn - this uses the package file dcolumn.sty to define two new types of column alignment for use in tables. If 
^ ] the usedcolumn option has been specified then, within a table, d{x} can be used to produce a 'flush right' decimal-aligned 
C/3 ■ column with x decimal places and . can be used to provide a decimal- aligned column centred on the decimal point (the 
, ^ , ' number of decimal places does not need to be specified for this column type). Note that the standard I^T^gX 'tools' packages 
, dcolumn.sty and array. sty are required in order to use the usedcolumn option. 

■ (iii) usenatbib - this uses Patrick Daly's natbib.sty package for cross-referencing. If the usenatbib option is specified, 
^ ' citations in the text should be in one of the following forms (or one of the additional forms documented within natbib . sty 

IT) '. itself). 
^ . 

. • \citet-[fce?/} produces text citations, e.g. Jones et al. (1990), 

■ • \citep{.keyy produces citations in parentheses, e.g. (Jones et al. 1990), 
• \citealt-[fce?/} produces citations with no parentheses, e.g. Jones et al. 1990. 



For three-author papers, a full author list can be forced by putting a * just before the {. To add notes within the citation, use 
the form \citep Ipre.reference.textl Lpost^reference^textlikey} (note that either of prc-reference^text and postjreference-text 
can be blank). 

Items in the reference list must be of the form 
\bibitem[\protect\citeauthoryear{ aMi/ior_names}-{j/ear}-] -Cfcej/} Text of reference ... 
for one-, two- and multi-author papers, or 

\bibitem [\protect\citeauthoryear{t/iree_aiithor_names3-{/trsLaiit/ior_etaZ}-{j/ear}-] {fcej/} Text of reference ... 
for three-author papers. 

Note that Patrick Daly's package natbib.sty is required in order to use the usenatbib option. 

We recommend that authors use natbib . sty as their standard cross-referencing package, because of the flexibility in citation 
style that it provides. 

(iv) usegraphicx - this loads the graphics package graphicx . sty, which authors can use to include figures in their papers. 
Note that the standard graphics package graphicx. sty is required in order to use the usegraphicx option. 

Other style files, or packages providing features such as graphics support, can be used in conjunction with mn.cls. To do 
this, the command \usepackage-[pacfca(;e_name} is used. 

An additional feature of the class file is that footnote symbols are now in the correct journal style (symbols for title page 
footnotes, superscript arable numbers for text footnotes). 

A BibTeX style file, mn2e.bst, is now available for authors who wish to use BiBTeX. It is recommended that this should 
be used in conjunction with natbib.sty. 

For general instructions on preparing Monthly Notices papers in ffl^rjX, please refer to the MNRAS M^^Xstyle guide for 
authors, available on the CTAN sites in subdirectory /tex-archive/macros/latex209/mnras. 



